!>\file micro_mg_utils.F90
!! This file contains process rates and utility functions used by the
!! MG microphysics.

!>\ingroup mg2mg3
!>\defgroup micro_mg_utils_mod Morrison-Gettelman MP utils Module
!! This module contains process rates and utility functions used by the MG 
!! microphysics.
!!
!! Original MG authors: Andrew Gettelman, Hugh Morrison
!! Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan
!!
!! Separated from MG 1.5 by B. Eaton.
!!
!! Separated module switched to MG 2.0 and further changes by S. Santos.
!!
!! Anning Cheng changed for FV3GFS 9/29/2017
!!              added ac_time as an input
!!
!! S. Moorthi - Feb 2018 : code optimization
!!
!! This version: https://svn-ccsm-models.cgd.ucar.edu/cam1/branch_tags/mg3_tags/mg3_33_cam5_4_153/
!!
!! for questions contact Hugh Morrison, Andrew Gettelman
!! e-mail: morrison@ucar.edu, andrew@ucar.edu
module micro_mg_utils

!--------------------------------------------------------------------------
!
! List of required external functions that must be supplied:
!   gamma --> standard mathematical gamma function (if gamma is an
!       intrinsic, define HAVE_GAMMA_INTRINSICS)
!
!--------------------------------------------------------------------------
!
! Constants that must be specified in the "init" method (module variables):
!
! kind            kind of reals (to verify correct linkage only) -
! gravit          acceleration due to gravity                    m s-2
! rair            dry air gas constant for air                   J kg-1 K-1
! rh2o            gas constant for water vapor                   J kg-1 K-1
! cpair           specific heat at constant pressure for dry air J kg-1 K-1
! tmelt           temperature of melting point for water         K
! latvap          latent heat of vaporization                    J kg-1
! latice          latent heat of fusion                          J kg-1
!
!--------------------------------------------------------------------------

! 8 byte real and integer
use machine,       only : r8 => kind_phys
use machine,       only : i8 => kind_phys
implicit none
private
save

public ::                           &
     micro_mg_utils_init,           &
     size_dist_param_liq,           &
     size_dist_param_basic,         &
     avg_diameter,                  &
     rising_factorial,              &
     ice_deposition_sublimation,    &
     sb2001v2_liq_autoconversion,   &
     sb2001v2_accre_cld_water_rain, &       
     kk2000_liq_autoconversion,     &
     ice_autoconversion,            &
     immersion_freezing,            &
     contact_freezing,              &
     snow_self_aggregation,         &
     accrete_cloud_water_snow,      &
     secondary_ice_production,      &
     accrete_rain_snow,             &
     heterogeneous_rain_freezing,   &
     accrete_cloud_water_rain,      &
     self_collection_rain,          &
     accrete_cloud_ice_snow,        &
     evaporate_sublimate_precip,    &
     bergeron_process_snow,         &
     liu_liq_autoconversion,        &
     gmao_ice_autoconversion,       &
     size_dist_param_ice,           &
!++ag
     graupel_collecting_snow,       &
     graupel_collecting_rain,       &
     graupel_collecting_cld_water,  &
     graupel_riming_liquid_snow,    &
     graupel_rain_riming_snow,      &
     graupel_rime_splintering,      &
     evaporate_sublimate_precip_graupel
!     graupel_sublimate_evap
!--ag


public :: MGHydrometeorProps

type :: MGHydrometeorProps
   ! Density (kg/m^3)
   real(r8) :: rho
   ! Information for size calculations.
   ! Basic calculation of mean size is:
   !     lambda = (shape_coef*nic/qic)^(1/eff_dim)
   ! Then lambda is constrained by bounds.
   real(r8) :: eff_dim
   real(r8) :: shape_coef
   real(r8) :: lambda_bounds(2)
   ! Minimum average particle mass (kg).
   ! Limit is applied at the beginning of the size distribution calculations.
   real(r8) :: min_mean_mass
end type MGHydrometeorProps

interface MGHydrometeorProps
   module procedure NewMGHydrometeorProps
end interface

type(MGHydrometeorProps), public :: mg_liq_props
type(MGHydrometeorProps), public :: mg_ice_props
type(MGHydrometeorProps), public :: mg_rain_props
type(MGHydrometeorProps), public :: mg_snow_props
!++ag
type(MGHydrometeorProps), public :: mg_graupel_props
!--ag

interface size_dist_param_liq
  module procedure size_dist_param_liq_vect
  module procedure size_dist_param_liq_line
end interface
interface size_dist_param_basic
  module procedure size_dist_param_basic_vect
  module procedure size_dist_param_basic_line
end interface

interface size_dist_param_ice
  module procedure size_dist_param_ice_vect
  module procedure size_dist_param_ice_line
end interface

!=================================================
! Public module parameters (mostly for MG itself)
!=================================================

!> Pi to 20 digits; more than enough to reach the limit of double precision.
real(r8), parameter, public :: pi = 3.14159265358979323846_r8

!> "One minus small number": number near unity for round-off issues.
!real(r8), parameter, public :: omsm   = 1._r8 - 1.e-5_r8
real(r8), parameter, public :: omsm   = 1._r8 - 1.e-6_r8

!> Smallest mixing ratio considered in microphysics.
real(r8), parameter, public :: qsmall = 1.e-18_r8

!> minimum allowed cloud fraction
 real(r8), parameter, public :: mincld = 0.000001_r8
!real(r8), parameter, public :: mincld = 0.0001_r8
!real(r8), parameter, public :: mincld = 0.0_r8

real(r8), parameter, public :: rhosn = 250._r8  !< bulk density snow
real(r8), parameter, public :: rhoi  = 500._r8  !< bulk density ice
real(r8), parameter, public :: rhow  = 1000._r8 !< bulk density liquid
real(r8), parameter, public :: rhows = 917._r8  !< bulk density water solid

!++ag
!Hail and Graupel (set in MG3)
real(r8), parameter, public :: rhog = 500._r8
real(r8), parameter, public :: rhoh = 400._r8
!--ag

! fall speed parameters, V = aD^b (V is in m/s)
! droplets
real(r8), parameter, public :: ac = 3.e7_r8
real(r8), parameter, public :: bc = 2._r8
! snow
real(r8), parameter, public :: as = 11.72_r8
real(r8), parameter, public :: bs = 0.41_r8
! cloud ice
real(r8), parameter, public :: ai = 700._r8
real(r8), parameter, public :: bi = 1._r8
! small cloud ice (r< 10 um) - sphere, bulk density
real(r8), parameter, public :: aj = ac*((rhoi/rhows)**(bc/3._r8))*rhows/rhow
real(r8), parameter, public :: bj = bc
! rain
real(r8), parameter, public :: ar = 841.99667_r8
real(r8), parameter, public :: br = 0.8_r8
!++ag
! graupel
real(r8), parameter, public :: ag = 19.3_r8
real(r8), parameter, public :: bg = 0.37_r8
! hail
real(r8), parameter, public :: ah = 114.5_r8
real(r8), parameter, public :: bh = 0.5_r8
!--ag

!> mass of new crystal due to aerosol freezing and growth (kg)
!! Make this consistent with the lower bound, to support UTLS and
!! stratospheric ice, and the smaller ice size limit.
real(r8), parameter, public :: mi0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3

!++ag
! mass of new graupel particle  (assume same as mi0 for now, may want to make bigger?)
!real(r8), parameter, public :: mg0 = 4._r8/3._r8*pi*rhoi*(1.e-6_r8)**3
!or set based on M2005:
real(r8), parameter, public :: mg0 = 1.6e-10_r8
! radius of contact nuclei
real(r8), parameter, public :: mmult = 4._r8/3._r8*pi*rhoi*(5.e-6_r8)**3
!--ag

!=================================================
! Private module parameters
!=================================================

! Signaling NaN bit pattern that represents a limiter that's turned off.
integer(i8), parameter :: limiter_off = int(Z'7FF1111111111111', i8)

! alternate threshold used for some in-cloud mmr
real(r8), parameter :: icsmall = 1.e-8_r8

! particle mass-diameter relationship
! currently we assume spherical particles for cloud ice/snow
! m = cD^d
! exponent
real(r8), parameter :: dsph = 3._r8

! Bounds for mean diameter for different constituents.
real(r8), parameter :: lam_bnd_rain(2) = 1._r8/[500.e-6_r8, 20.e-6_r8]
real(r8), parameter :: lam_bnd_snow(2) = 1._r8/[2000.e-6_r8, 10.e-6_r8]

! Minimum average mass of particles.
real(r8), parameter :: min_mean_mass_liq = 1.e-20_r8
real(r8), parameter :: min_mean_mass_ice = 1.e-20_r8

! ventilation parameters
! for snow
real(r8), parameter :: f1s = 0.86_r8
real(r8), parameter :: f2s = 0.28_r8
! for rain
real(r8), parameter :: f1r = 0.78_r8
real(r8), parameter :: f2r = 0.308_r8

! collection efficiencies
! aggregation of cloud ice and snow
!real(r8), parameter :: eii = 0.5_r8
!real(r8), parameter :: eii = 0.1_r8
 real(r8), parameter :: eii = 0.2_r8
!++ag
! collection efficiency, ice-droplet collisions
real(r8), parameter, public :: ecid = 0.7_r8
! collection efficiency between droplets/rain and snow/rain
real(r8), parameter, public :: ecr = 1.0_r8
!--ag

! immersion freezing parameters, bigg 1953
real(r8), parameter :: bimm = 100._r8
real(r8), parameter :: aimm = 0.66_r8

! Mass of each raindrop created from autoconversion.
real(r8), parameter :: droplet_mass_25um = 4._r8/3._r8*pi*rhow*(25.e-6_r8)**3
real(r8), parameter :: droplet_mass_40um = 4._r8/3._r8*pi*rhow*(40.e-6_r8)**3, &
                       droplet_mass_40umi = 1._r8/droplet_mass_40um

!=========================================================
! Constants set in initialization
!=========================================================

! Set using arguments to micro_mg_init
real(r8) :: rv          ! water vapor gas constant
real(r8) :: cpp         ! specific heat of dry air
real(r8) :: tmelt       ! freezing point of water (K)

real(r8) :: ra        ! dry air gas constant

! latent heats of:
real(r8) :: xxlv        ! vaporization
real(r8) :: xlf         ! freezing
real(r8) :: xxls        ! sublimation

! additional constants to help speed up code
real(r8) :: gamma_bs_plus3
real(r8) :: gamma_half_br_plus5
real(r8) :: gamma_half_bs_plus5
!++ag
real(r8) :: gamma_2bs_plus2
!--ag
!
real(r8), parameter :: zero = 0._r8, one = 1._r8,  two = 2._r8,  three = 3._r8, &
                       four = 4._r8, five = 5._r8, six  = 6._r8, pio6 = pi/six, &
                       pio3 = pi/three, half = 0.5_r8, oneo3 = one/three,       &
                       twopi = pi + pi

!=========================================================
! Utilities that are cheaper if the compiler knows that
! some argument is an integer.
!=========================================================

!>\ingroup micro_mg_utils_mod
interface rising_factorial
   module procedure rising_factorial_r8
   module procedure rising_factorial_integer
end interface rising_factorial

!>\ingroup micro_mg_utils_mod
interface var_coef
   module procedure var_coef_r8
   module procedure var_coef_integer
end interface var_coef

!==========================================================================
contains
!==========================================================================

!>\ingroup micro_mg_utils_mod 
!! Initialize module variables.
!
! "kind" serves no purpose here except to check for unlikely linking
! issues; always pass in the kind for a double precision real.
!
!
! Check the list at the top of this module for descriptions of all other
! arguments.
subroutine micro_mg_utils_init( kind, rair, rh2o, cpair, tmelt_in, latvap, &
                                latice, dcs)
!                               latice, dcs, errstring)

  integer,  intent(in)  :: kind
!++ag
  real(r8), intent(in)  :: rair
!--ag
  real(r8), intent(in)  :: rh2o
  real(r8), intent(in)  :: cpair
  real(r8), intent(in)  :: tmelt_in
  real(r8), intent(in)  :: latvap
  real(r8), intent(in)  :: latice
  real(r8), intent(in)  :: dcs


  ! Name this array to workaround an XLF bug (otherwise could just use the
  ! expression that sets it).
  real(r8) :: ice_lambda_bounds(2)

  !-----------------------------------------------------------------------


  ! declarations for MG code (transforms variable names)

  rv    = rh2o          ! water vapor gas constant
  cpp   = cpair         ! specific heat of dry air
  tmelt = tmelt_in
  !++ag
  ra = rair             ! dry air gas constant
  !--ag

  ! latent heats

  xxlv = latvap         ! latent heat vaporization
  xlf  = latice         ! latent heat freezing
  xxls = xxlv + xlf     ! latent heat of sublimation

  ! Define constants to help speed up code (this limits calls to gamma function)
  gamma_bs_plus3      = gamma(three+bs)
  gamma_half_br_plus5 = gamma((five+br)*half)
  gamma_half_bs_plus5 = gamma((five+bs)*half)
!++ag
  gamma_2bs_plus2     = gamma(bs+bs+two)
!--ag


  ! Don't specify lambda bounds for cloud liquid, as they are determined by
  ! pgam dynamically.
  mg_liq_props = MGHydrometeorProps(rhow, dsph, min_mean_mass=min_mean_mass_liq)

  ! Mean ice diameter can not grow bigger than twice the autoconversion
  ! threshold for snow.
  ice_lambda_bounds = one/[two*dcs, 1.e-6_r8]

  mg_ice_props      = MGHydrometeorProps(rhoi, dsph, &
                                        ice_lambda_bounds, min_mean_mass_ice)

  mg_rain_props     = MGHydrometeorProps(rhow,  dsph, lam_bnd_rain)
  mg_snow_props     = MGHydrometeorProps(rhosn, dsph, lam_bnd_snow)
!++ag
  mg_graupel_props  = MGHydrometeorProps(rhog, dsph, lam_bnd_snow)
!--ag

end subroutine micro_mg_utils_init

!>\ingroup micro_mg_utils_mod
!! Constructor for a constituent property object.
function NewMGHydrometeorProps(rho, eff_dim, lambda_bounds, min_mean_mass) &
     result(res)
  real(r8), intent(in) :: rho, eff_dim
  real(r8), intent(in), optional :: lambda_bounds(2), min_mean_mass
  type(MGHydrometeorProps) :: res

  res%rho = rho
  res%eff_dim = eff_dim
  if (present(lambda_bounds)) then
     res%lambda_bounds = lambda_bounds
  else
     res%lambda_bounds = no_limiter()
  end if
  if (present(min_mean_mass)) then
     res%min_mean_mass = min_mean_mass
  else
     res%min_mean_mass = no_limiter()
  end if

  res%shape_coef = rho * pio6 * gamma(eff_dim+one)

end function NewMGHydrometeorProps

!========================================================================
!FORMULAS
!========================================================================

! Use gamma function to implement rising factorial extended to the reals.
pure function rising_factorial_r8(x, n) result(res)
  real(r8), intent(in) :: x, n
  real(r8) :: res

  res = gamma(x+n) / gamma(x)

end function rising_factorial_r8

! Rising factorial can be performed much cheaper if n is a small integer.
pure function rising_factorial_integer(x, n) result(res)
  real(r8), intent(in) :: x
  integer,  intent(in) :: n
  real(r8) :: res

  integer  :: i
  real(r8) :: factor

  res    = one
  factor = x

  do i = 1, n
     res    = res * factor
     factor = factor + one
  end do

end function rising_factorial_integer

! Calculate correction due to latent heat for evaporation/sublimation
elemental function calc_ab(t, qv, xxl) result(ab)
  real(r8), intent(in) :: t     ! Temperature
  real(r8), intent(in) :: qv    ! Saturation vapor pressure
  real(r8), intent(in) :: xxl   ! Latent heat

  real(r8) :: ab

  real(r8) :: dqsdt

  dqsdt = xxl*qv / (rv*t*t)
  ab = one + dqsdt*xxl/cpp

end function calc_ab

!>\ingroup micro_mg_utils_mod
!! get cloud droplet size distribution parameters
elemental subroutine size_dist_param_liq_line(props, qcic, ncic, rho, pgam, lamc)
  type(MGHydrometeorProps), intent(in) :: props
  real(r8), intent(in)     :: qcic
  real(r8), intent(inout)  :: ncic
  real(r8), intent(in)     :: rho

  real(r8), intent(out)    :: pgam
  real(r8), intent(out)    :: lamc
  real(r8)                 :: xs

  type(MGHydrometeorProps) :: props_loc
  logical, parameter       :: liq_gmao=.true.

  if (qcic > qsmall) then

     ! Local copy of properties that can be modified.
     ! (Elemental routines that operate on arrays can't modify scalar
     ! arguments.)
     props_loc = props

     ! Get pgam from fit to observations of martin et al. 1994

     if (liq_gmao) then
       pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8
       ! Anning modified lamc
       if ((ncic > 1.0e-3_r8) .and. (qcic > 1.0e-11_r8)) then
         xs = 0.07_r8*(1000._r8*qcic/ncic) ** (-0.14_r8)
       else
         xs = 1.2_r8
       end if

        xs = max(min(xs, 1.7_r8), 1.1_r8)
        xs = xs*xs*xs
        xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8
        pgam = sqrt(xs)
     else

       pgam = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic*rho)
!      pgam = 0.0005714_r8*1.e-6_r8*ncic*rho + 0.2714_r8
     endif
     pgam = one / (pgam*pgam) - one
     pgam = max(pgam, two)

     ! Set coefficient for use in size_dist_param_basic.
     ! The 3D case is so common and optimizable that we specialize it:
     if (props_loc%eff_dim == three) then
        props_loc%shape_coef = pio6 * props_loc%rho * &
                               rising_factorial(pgam+one, 3)
     else
        props_loc%shape_coef = pio6 * props_loc%rho * &
                               rising_factorial(pgam+one, props_loc%eff_dim)
     end if

     ! Limit to between 2 and 50 microns mean size.
     props_loc%lambda_bounds = (pgam+one) * one/[50.e-6_r8, 2.e-6_r8]

     call size_dist_param_basic(props_loc, qcic, ncic, lamc)

  else
     ! pgam not calculated in this case, so set it to a value likely to
     ! cause an error if it is accidentally used
     ! (gamma function undefined for negative integers)
     pgam = -100._r8
     lamc = zero
  end if

end subroutine size_dist_param_liq_line

!>\ingroup micro_mg_utils_mod
!! This subroutine gets cloud droplet size distribution parameters
subroutine size_dist_param_liq_vect(props, qcic, ncic, rho, pgam, lamc, mgncol)

  type(mghydrometeorprops),    intent(in)    :: props
  integer,                     intent(in)    :: mgncol
  real(r8), dimension(mgncol), intent(in)    :: qcic
  real(r8), dimension(mgncol), intent(inout) :: ncic
  real(r8), dimension(mgncol), intent(in)    :: rho
  real(r8), dimension(mgncol), intent(out)   :: pgam
  real(r8), dimension(mgncol), intent(out)   :: lamc
  real(r8)                                   :: xs
  type(mghydrometeorprops)                   :: props_loc
  logical, parameter                         :: liq_gmao=.true.
  integer :: i

  do i=1,mgncol
     if (qcic(i) > qsmall) then
        ! Local copy of properties that can be modified.
        ! (Elemental routines that operate on arrays can't modify scalar
        ! arguments.)
        props_loc = props
        ! Get pgam from fit to observations of martin et al. 1994

        if (liq_gmao) then
          pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8
          if ((ncic(i) > 1.0e-3_r8) .and. (qcic(i) > 1.0e-11_r8)) then
            xs = 0.07_r8*(1000._r8*qcic(i)/ncic(i)) **(-0.14_r8)
          else
            xs = 1.2_r8
          end if

          xs = max(min(xs, 1.7_r8), 1.1_r8)
          xs = xs*xs*xs
          xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) - 4.0_r8)/8.0_r8
          pgam(i) = sqrt(xs)
        else
          pgam(i) = one - 0.7_r8 * exp(-0.008_r8*1.e-6_r8*ncic(i)*rho(i))
!         pgam(i) = 0.0005714_r8*1.e-6_r8*ncic(i)*rho(i) + 0.2714_r8
        endif

        pgam(i) = one/(pgam(i)*pgam(i)) - one
        pgam(i) = max(pgam(i), two)
     endif
  enddo
  do i=1,mgncol
     if (qcic(i) > qsmall) then
        ! Set coefficient for use in size_dist_param_basic.
        ! The 3D case is so common and optimizable that we specialize
        ! it:
        if (props_loc%eff_dim == three) then
           props_loc%shape_coef = pio6 * props_loc%rho * &
                rising_factorial(pgam(i)+one, 3)
        else
           props_loc%shape_coef = pio6 * props_loc%rho * &
                rising_factorial(pgam(i)+one, props_loc%eff_dim)
        end if
        ! Limit to between 2 and 50 microns mean size.
        props_loc%lambda_bounds(1) = (pgam(i)+one) / 50.e-6_r8
        props_loc%lambda_bounds(2) = (pgam(i)+one) / 2.e-6_r8
        call size_dist_param_basic(props_loc, qcic(i), ncic(i), lamc(i))
     endif
  enddo
  do i=1,mgncol
     if (qcic(i) <= qsmall) then
        ! pgam not calculated in this case, so set it to a value likely to
        ! cause an error if it is accidentally used
        ! (gamma function undefined for negative integers)
        pgam(i) = -100._r8
        lamc(i) = zero
     end if
  enddo

end subroutine size_dist_param_liq_vect

!>\ingroup micro_mg_utils_mod
!! Basic routine for getting size distribution parameters.
elemental subroutine size_dist_param_basic_line(props, qic, nic, lam, n0)
  type(MGHydrometeorProps), intent(in) :: props
  real(r8), intent(in)    :: qic
  real(r8), intent(inout) :: nic

  real(r8), intent(out)   :: lam
  real(r8), intent(out), optional :: n0

  if (qic > qsmall) then

     ! add upper limit to in-cloud number concentration to prevent
     ! numerical error
     if (limiter_is_on(props%min_mean_mass)) then
        nic = min(nic, qic / props%min_mean_mass)
     end if

     ! lambda = (c n/q)^(1/d)
     lam = (props%shape_coef * nic/qic)**(one/props%eff_dim)

     ! check for slope
     ! adjust vars
     if (lam < props%lambda_bounds(1)) then
        lam = props%lambda_bounds(1)
        nic = lam**(props%eff_dim) * qic/props%shape_coef
     else if (lam > props%lambda_bounds(2)) then
        lam = props%lambda_bounds(2)
        nic = lam**(props%eff_dim) * qic/props%shape_coef
     end if

  else
     lam = zero
  end if

  if (present(n0)) n0 = nic * lam

end subroutine size_dist_param_basic_line

!>\ingroup micro_mg_utils_mod
!! This subroutine calculates
subroutine size_dist_param_basic_vect(props, qic, nic, lam, mgncol, n0)

  type (mghydrometeorprops),   intent(in)    :: props
  integer,                     intent(in)    :: mgncol
  real(r8), dimension(mgncol), intent(in)    :: qic
  real(r8), dimension(mgncol), intent(inout) :: nic
  real(r8), dimension(mgncol), intent(out)   :: lam
  real(r8), dimension(mgncol), intent(out), optional :: n0
  integer :: i
  do i=1,mgncol

     if (qic(i) > qsmall) then

        ! add upper limit to in-cloud number concentration to prevent
        ! numerical error
        if (limiter_is_on(props%min_mean_mass)) then
           nic(i) = min(nic(i), qic(i) / props%min_mean_mass)
        end if

        ! lambda = (c n/q)^(1/d)
        lam(i) = (props%shape_coef * nic(i)/qic(i))**(one/props%eff_dim)

        ! check for slope
        ! adjust vars
        if (lam(i) < props%lambda_bounds(1)) then
           lam(i) = props%lambda_bounds(1)
           nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
        else if (lam(i) > props%lambda_bounds(2)) then
           lam(i) = props%lambda_bounds(2)
           nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
        end if

     else
        lam(i) = zero
     end if

  enddo

  if (present(n0)) n0 = nic * lam

end subroutine size_dist_param_basic_vect

!>\ingroup micro_mg_utils_mod
!! ice routine for getting size distribution parameters.
elemental subroutine size_dist_param_ice_line(props, qic, nic, lam, n0)
  type(MGHydrometeorProps), intent(in) :: props
  real(r8), intent(in) :: qic
  real(r8), intent(inout) :: nic

  real(r8), intent(out) :: lam
  real(r8):: miu_ice,tx1,tx2, aux
  real(r8), intent(out), optional :: n0
  logical, parameter  ::  ice_sep=.true.

  if (qic > qsmall) then

     ! add upper limit to in-cloud number concentration to prevent
     ! numerical error
     if (limiter_is_on(props%min_mean_mass)) then
        nic = min(nic, qic / props%min_mean_mass)
     end if

     ! lambda = (c n/q)^(1/d)
     lam = (props%shape_coef * nic/qic)**(1._r8/props%eff_dim)
     if (ice_sep) then
       miu_ice = max(min(0.008_r8*(lam*0.01)**0.87_r8, 10.0_r8), 0.1_r8)
       tx1 = 1.0_r8 + miu_ice
       tx2 = 1.0_r8 / gamma(tx1)
       aux = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8)
       lam = lam*aux
     else
       aux = 1.0_r8
       tx1 = 1.0_r8
       tx2 = 1.0_r8
     end if
     if (present(n0)) n0 = nic * lam**tx1*tx2

     ! check for slope
     ! adjust vars
     if (lam < props%lambda_bounds(1)*aux) then
        lam = props%lambda_bounds(1)
        nic = lam**(props%eff_dim) * qic/props%shape_coef
        if (present(n0)) n0 = nic * lam
     else if (lam > props%lambda_bounds(2)*aux) then
        lam = props%lambda_bounds(2)
        nic = lam**(props%eff_dim) * qic/props%shape_coef
        if (present(n0)) n0 = nic * lam
     end if

  else
     lam = 0.0_r8
  end if


end subroutine size_dist_param_ice_line

!>\ingroup micro_mg_utils_mod
!! This subroutine
subroutine size_dist_param_ice_vect(props, qic, nic, lam, mgncol, n0)

  type (mghydrometeorprops), intent(in) :: props
  integer,                          intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: qic
  real(r8), dimension(mgncol), intent(inout) :: nic
  real(r8), dimension(mgncol), intent(out) :: lam
  real(r8), dimension(mgncol), intent(out), optional :: n0
  real(r8) :: miu_ice,tx1,tx2, aux
  integer :: i
  logical, parameter  ::  ice_sep=.true.
  do i=1,mgncol

     if (qic(i) > qsmall) then

        ! add upper limit to in-cloud number concentration to prevent
        ! numerical error
        if (limiter_is_on(props%min_mean_mass)) then
           nic(i) = min(nic(i), qic(i) / props%min_mean_mass)
        end if

        ! lambda = (c n/q)^(1/d)
        lam(i) = (props%shape_coef * nic(i)/qic(i))**(1._r8/props%eff_dim)
        if (ice_sep) then
           miu_ice = max(min(0.008_r8*(lam(i)*0.01)**0.87_r8, 10.0_r8), 0.1_r8)
           tx1    = 1.0_r8 + miu_ice
           tx2    = 1.0_r8 / gamma(tx1)
           aux    = (gamma(tx1+3.0_r8)*tx2) ** (1.0_r8/3.0_r8)
           lam(i) = lam(i)*aux
        else
           aux = 1.0_r8
           tx1 = 1.0_r8
           tx2 = 1.0_r8
        end if
        if (present(n0)) n0(i) = nic(i) * lam(i)**tx1*tx2

        ! check for slope
        ! adjust vars
        if (lam(i) < props%lambda_bounds(1)*aux) then
           lam(i) = props%lambda_bounds(1)
           nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
           if (present(n0)) n0(i) = nic(i) * lam(i)
        else if (lam(i) > props%lambda_bounds(2)*aux) then
           lam(i) = props%lambda_bounds(2)
           nic(i) = lam(i)**(props%eff_dim) * qic(i)/props%shape_coef
           if (present(n0)) n0(i) = nic(i) * lam(i)
        end if

     else
        lam(i) = 0.0_r8
     end if

  enddo

end subroutine size_dist_param_ice_vect

!>\ingroup micro_mg_utils_mod
!> Finds the average diameter of particles given their density, and
!! mass/number concentrations in the air.
!! Assumes that diameter follows an exponential distribution.
real(r8) elemental function avg_diameter(q, n, rho_air, rho_sub)
  real(r8), intent(in) :: q         !< mass mixing ratio
  real(r8), intent(in) :: n         !< number concentration (per volume)
  real(r8), intent(in) :: rho_air   !< local density of the air
  real(r8), intent(in) :: rho_sub   !< density of the particle substance

  avg_diameter = (pi * rho_sub * n/(q*rho_air))**(-oneo3)

end function avg_diameter

!>\ingroup mg2mg3
!> Finds a coefficient for process rates based on the relative variance
!! of cloud water.
elemental function var_coef_r8(relvar, a) result(res)
  real(r8), intent(in) :: relvar
  real(r8), intent(in) :: a
  real(r8) :: res

  res = rising_factorial(relvar, a) / relvar**a

end function var_coef_r8

!>\ingroup mg2mg3
!> Finds a coefficient for process rates based on the relative variance
!! of cloud water.
elemental function var_coef_integer(relvar, a) result(res)
  real(r8), intent(in) :: relvar
  integer, intent(in) :: a
  real(r8) :: res

  res = rising_factorial(relvar, a) / relvar**a

end function var_coef_integer

!========================================================================
!MICROPHYSICAL PROCESS CALCULATIONS
!========================================================================
!========================================================================
!>\ingroup micro_mg_utils_mod
!! Initial ice deposition and sublimation loop.
!! Run before the main loop
!! This subroutine written by Peter Caldwell
subroutine ice_deposition_sublimation(t, qv, qi, ni,           &
                                      icldm, rho, dv,qvl, qvi, &
                                      berg, vap_dep, ice_sublim, mgncol)

  !INPUT VARS:
  !===============================================
  integer,  intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t
  real(r8), dimension(mgncol), intent(in) :: qv
  real(r8), dimension(mgncol), intent(in) :: qi
  real(r8), dimension(mgncol), intent(in) :: ni
  real(r8), dimension(mgncol), intent(in) :: icldm
  real(r8), dimension(mgncol), intent(in) :: rho
  real(r8), dimension(mgncol), intent(in) :: dv
  real(r8), dimension(mgncol), intent(in) :: qvl
  real(r8), dimension(mgncol), intent(in) :: qvi

  !OUTPUT VARS:
  !===============================================
  real(r8), dimension(mgncol), intent(out) :: vap_dep !ice deposition (cell-ave value)
  real(r8), dimension(mgncol), intent(out) :: ice_sublim !ice sublimation (cell-ave value)
  real(r8), dimension(mgncol), intent(out) :: berg !bergeron enhancement (cell-ave value)

  !INTERNAL VARS:
  !===============================================
  real(r8) :: ab
  real(r8) :: epsi
  real(r8) :: qiic
  real(r8) :: niic
  real(r8) :: lami
  real(r8) :: n0i
  real(r8) :: tx1
  integer  :: i

  do i=1,mgncol
     if (qi(i)>=qsmall) then

        !GET IN-CLOUD qi, ni
        !===============================================
        tx1  = one / icldm(i)
        qiic = qi(i) * tx1
        niic = ni(i) * tx1

        !Compute linearized condensational heating correction
        ab = calc_ab(t(i), qvi(i), xxls)
        !Get slope and intercept of gamma distn for ice.
!       call size_dist_param_basic(mg_ice_props, qiic, niic, lami, n0i)
        call size_dist_param_ice(mg_ice_props, qiic, niic, lami, n0i)
        !Get depletion timescale=1/eps
        epsi = twopi*n0i*rho(i)*Dv(i)/(lami*lami)

        !Compute deposition/sublimation
        vap_dep(i) = epsi/ab*(qv(i) - qvi(i))

        !Make this a grid-averaged quantity
        vap_dep(i) = vap_dep(i)*icldm(i)

        !Split into deposition or sublimation.
        if (t(i) < tmelt .and. vap_dep(i) > zero) then
           ice_sublim(i) = zero
        else
        ! make ice_sublim negative for consistency with other evap/sub processes
           ice_sublim(i) = min(vap_dep(i), zero)
           vap_dep(i) = zero
        end if

        !sublimation occurs @ any T. Not so for berg.
        if (t(i) < tmelt) then

           !Compute bergeron rate assuming cloud for whole step.
           berg(i) = max(epsi/ab*(qvl(i) - qvi(i)), zero)
        else !T>frz
           berg(i) = zero
        end if !T<frz

     else !where qi<qsmall
        berg(i)       = zero
        vap_dep(i)    = zero
        ice_sublim(i) = zero
     end if !qi>qsmall
  enddo
end subroutine ice_deposition_sublimation

!========================================================================
!>\ingroup micro_mg_utils_mod
!! autoconversion of cloud liquid water to rain
!! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
!! minimum qc of 1 x 10^-8 prevents floating point error
subroutine kk2000_liq_autoconversion(microp_uniform, qcic, &
                            ncic, rho, relvar, prc, nprc, nprc1, mgncol)

  integer, intent(in) :: mgncol
  logical, intent(in) :: microp_uniform

  real(r8), dimension(mgncol), intent(in) :: qcic
  real(r8), dimension(mgncol), intent(in) :: ncic
  real(r8), dimension(mgncol), intent(in) :: rho

  real(r8), dimension(mgncol), intent(in) :: relvar

  real(r8), dimension(mgncol), intent(out) :: prc
  real(r8), dimension(mgncol), intent(out) :: nprc
  real(r8), dimension(mgncol), intent(out) :: nprc1

  real(r8), dimension(mgncol) :: prc_coef
  integer :: i

  ! Take variance into account, or use uniform value.
  if (.not. microp_uniform) then
     prc_coef(:) = var_coef(relvar(:), 2.47_r8)
  else
     prc_coef(:) = one
  end if

  do i=1,mgncol
     if (qcic(i) >= icsmall) then

        ! nprc is increase in rain number conc due to autoconversion
        ! nprc1 is decrease in cloud droplet conc due to autoconversion

        ! assume exponential sub-grid distribution of qc, resulting in additional
        ! factor related to qcvar below
        ! switch for sub-columns, don't include sub-grid qc

        prc(i) = prc_coef(i) * &
             0.01_r8 * 1350._r8 * qcic(i)**2.47_r8 * (ncic(i)*1.e-6_r8*rho(i))**(-1.1_r8)
        nprc(i)  = prc(i) * (one/droplet_mass_25um)
        nprc1(i) = prc(i)*ncic(i)/qcic(i)

     else
        prc(i)   = zero
        nprc(i)  = zero
        nprc1(i) = zero
     end if
  enddo
end subroutine kk2000_liq_autoconversion
  
  !========================================================================
!>\ingroup micro_mg_utils_mod
!! This subroutine
subroutine sb2001v2_liq_autoconversion(pgam,qc,nc,qr,rho,relvar,au,nprc,nprc1,mgncol)
  !
  ! ---------------------------------------------------------------------
  ! AUTO_SB:  calculates the evolution of mass- and number mxg-ratio for
  ! drizzle drops due to autoconversion. The autoconversion rate assumes
  ! f(x)=A*x**(nu_c)*exp(-Bx) in drop MASS x. 

  ! Code from Hugh Morrison, Sept 2014

  ! autoconversion
  ! use simple lookup table of dnu values to get mass spectral shape parameter
  ! equivalent to the size spectral shape parameter pgam
    
  integer, intent(in) :: mgncol  
  
  real(r8), dimension(mgncol), intent (in)    :: pgam
  real(r8), dimension(mgncol), intent (in)    :: qc    ! = qc (cld water mixing ratio)
  real(r8), dimension(mgncol), intent (in)    :: nc    ! = nc (cld water number conc /kg)    
  real(r8), dimension(mgncol), intent (in)    :: qr    ! = qr (rain water mixing ratio)
  real(r8), dimension(mgncol), intent (in)    :: rho   ! = rho : density profile
  real(r8), dimension(mgncol), intent (in)    :: relvar 
  
  real(r8), dimension(mgncol), intent (out)   :: au    ! = prc autoconversion rate
  real(r8), dimension(mgncol), intent (out)   :: nprc1 ! = number tendency
  real(r8), dimension(mgncol), intent (out)   :: nprc  ! = number tendency fixed size for rain
 
  ! parameters for droplet mass spectral shape, 
  ! used by Seifert and Beheng (2001)                             
  ! warm rain scheme only (iparam = 1)                                                                        
  real(r8), parameter :: dnu(16) = [0._r8,-0.557_r8,-0.430_r8,-0.307_r8, & 
     -0.186_r8,-0.067_r8,0.050_r8,0.167_r8,0.282_r8,0.397_r8,0.512_r8,   &
     0.626_r8,0.739_r8,0.853_r8,0.966_r8,0.966_r8]

  ! parameters for Seifert and Beheng (2001) autoconversion/accretion                                         
  real(r8), parameter :: kc  = 9.44e9_r8
  real(r8), parameter :: kr  = 5.78e3_r8
  real(r8), parameter :: auf = kc / (20._r8*2.6e-7_r8) * 1000._r8, &
                         con_nprc1 = two/2.6e-7_r8*1000._r8
  real(r8) :: dum, dum1, nu, pra_coef, tx1, tx2, tx3, tx4
  integer  :: dumi, i

  do i=1,mgncol

    pra_coef = var_coef(relvar(i), 2.47_r8)
     if (qc(i) > qsmall) then
       dumi = max(1, min(int(pgam(i)), 15))
       nu   = dnu(dumi) + (dnu(dumi+1)-dnu(dumi))* (pgam(i)-dumi)

       !Anning fixed a bug here for FV3GFS 10/13/2017
       dum  = max(one-qc(i)/(qc(i)+qr(i)), zero)  
       tx1  = dum**0.68_r8
       tx2  = one - tx1
       dum1 = 600._r8 * tx1 * tx2 * tx2 * tx2      ! Moorthi
!      dum1 = 600._r8*dum**0.68_r8*(one-dum**0.68_r8)**3

       tx1 = nu + one
       tx2 = 0.001_r8 * rho(i) * qc(i)
       tx3 = tx2 * tx2 / (rho(i)*nc(i)*1.e-6_r8)
       tx2 = tx3 * tx3
       tx3 = one - dum
       au(i) = auf * (nu+two) * (nu+four) * tx2       &
             * (one+dum1/(tx3*tx3)) / (tx1*tx1*rho(i))

!      au(i) = kc/(20._r8*2.6e-7_r8)* &
!        (nu+2._r8)*(nu+4._r8)/(nu+1._r8)**2._r8* &
!        (rho(i)*qc(i)/1000._r8)**4._r8/(rho(i)*nc(i)/1.e6_r8)**2._r8* &
!        (1._r8+dum1/(1._r8-dum)**2)*1000._r8 / rho(i)

!      nprc1(i) = au(i) * two / 2.6e-7_r8 * 1000._r8
!      nprc(i)  = au(i) / droplet_mass_40um
       nprc1(i) = au(i) * con_nprc1
       nprc(i)  = au(i) * droplet_mass_40umi
     else
       au(i)    = zero
       nprc1(i) = zero
       nprc(i)  = zero
     end if
  
  enddo

  end subroutine sb2001v2_liq_autoconversion 
  
!========================================================================
!>\ingroup micro_mg_utils_mod
!!  Anning Cheng 10/5/2017 add Liu et al. autoconversion
  subroutine liu_liq_autoconversion(pgam,qc,nc,qr,rho,relvar,  &
                                    au,nprc,nprc1,mgncol)



       integer, intent(in) :: mgncol

       real(r8), dimension(mgncol), intent (in) :: pgam
       real(r8), dimension(mgncol), intent (in) :: qc
       real(r8), dimension(mgncol), intent (in) :: nc
       real(r8), dimension(mgncol), intent (in) :: qr
       real(r8), dimension(mgncol), intent (in) :: rho
       real(r8), dimension(mgncol), intent (in) :: relvar

       real(r8), dimension(mgncol), intent (out) :: au
       real(r8), dimension(mgncol), intent (out) :: nprc1
       real(r8), dimension(mgncol), intent (out) :: nprc
       real(r8) :: xs,lw, nw, beta6
!      real(r8), parameter :: dcrit=1.0e-6, miu_disp=1.
!      real(r8), parameter :: dcrit=1.0e-3, miu_disp=1.
       real(r8), parameter :: dcrit = 2.0e-3, miu_disp = 0.8,     &
                              con_nprc1 = two/2.6e-7_r8*1000._r8
       integer ::  i

       do i=1,mgncol
         if (qc(i) > qsmall) then
           xs    = one / (one+pgam(i))
           beta6 = (one+three*xs)*(one+four*xs)*(one+five*xs)  &
                 / ((one+xs)*(one+xs+xs))
           LW    = 1.0e-3_r8 * qc(i) * rho(i)
           NW    = nc(i) * rho(i) * 1.e-6_r8

           xs    = min(20.0_r8, 1.03e16_r8*(LW*LW)/(NW*SQRT(NW)))
           au(i) = 1.1e10_r8*beta6*LW*LW*LW                &
                 * (one-exp(-(xs**miu_disp))) / NW
           au(i) = au(i)*1.0e3_r8/rho(i)
           au(i) = au(i) * gamma(two+relvar(i))          &
                 / (gamma(relvar(i))*(relvar(i)*relvar(i)))

           au(i)   = au(i) * dcrit
!          nprc1(i)= au(i) * (two/2.6e-7_r8*1000._r8)
!          nprc(i) = au(i) / droplet_mass_40um
           nprc1(i)= au(i) * con_nprc1
           nprc(i) = au(i) * droplet_mass_40umi
         else
           au(i)    = zero
           nprc1(i) = zero
           nprc(i)  = zero
         end if
       enddo

  end subroutine liu_liq_autoconversion


!========================================================================
!SB2001 Accretion V2
!>\ingroup micro_mg_utils_mod
subroutine sb2001v2_accre_cld_water_rain(qc,nc,qr,rho,relvar,pra,npra,mgncol)
  !
  ! ---------------------------------------------------------------------
  ! ACCR_SB calculates the evolution of mass mxng-ratio due to accretion
  ! and self collection following Seifert & Beheng (2001).  
  !
  
  integer, intent(in) :: mgncol
  
  real(r8), dimension(mgncol), intent (in)    :: qc  ! = qc (cld water mixing ratio)
  real(r8), dimension(mgncol), intent (in)    :: nc  ! = nc (cld water number conc /kg)    
  real(r8), dimension(mgncol), intent (in)    :: qr  ! = qr (rain water mixing ratio)
  real(r8), dimension(mgncol), intent (in)    :: rho ! = rho : density profile
  real(r8), dimension(mgncol), intent (in)    :: relvar

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: pra  ! MMR
  real(r8), dimension(mgncol), intent(out) :: npra ! Number

  ! parameters for Seifert and Beheng (2001) autoconversion/accretion                                         
  real(r8), parameter :: kc = 9.44e9_r8
  real(r8), parameter :: kr = 5.78e3_r8

  real(r8) :: dum, dum1, tx1, tx2
  integer :: i

  ! accretion

  do i =1,mgncol

    if (qc(i) > qsmall) then
      dum     = one - qc(i)/(qc(i)+qr(i))
      tx1     = dum / (dum+5.e-4_r8)
      dum1    = tx1 * tx1
      dum1    = dum1 * dum1
      pra(i)  = kr*rho(i)*0.001_r8*qc(i)*qr(i)*dum1

      npra(i) = pra(i) * nc(i) / qc(i)

!     npra(i) = pra(i)*rho(i)*0.001_r8*(nc(i)*rho(i)*1.e-6_r8)/ &
!              (qc(i)*rho(i)*0.001_r8)*1.e6_r8 / rho(i)
    else
      pra(i)  = zero
      npra(i) = zero
    end if 
  
  enddo
 
  end subroutine sb2001v2_accre_cld_water_rain   

!========================================================================
! Autoconversion of cloud ice to snow
! similar to Ferrier (1994)
!>\ingroup micro_mg_utils_mod
!! Autoconversion of cloud ice to snow
!! similar to Ferrier (1994)
subroutine ice_autoconversion(t, qiic, lami, n0i, dcs, ac_time, prci, nprci, mgncol)

  integer, intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t
  real(r8), dimension(mgncol), intent(in) :: qiic
  real(r8), dimension(mgncol), intent(in) :: lami
  real(r8), dimension(mgncol), intent(in) :: n0i
  real(r8),                    intent(in) :: dcs
  real(r8), dimension(mgncol), intent(in) :: ac_time 

  real(r8), dimension(mgncol), intent(out) :: prci
  real(r8), dimension(mgncol), intent(out) :: nprci

  ! Assume autoconversion timescale of 180 seconds.

  ! Average mass of an ice particle.
  real(r8) :: m_ip
  ! Ratio of autoconversion diameter to average diameter.
  real(r8) :: d_rat
  integer :: i

  do i=1,mgncol
     if (t(i) <= tmelt .and. qiic(i) >= qsmall) then

        d_rat = lami(i)*dcs

        ! Rate of ice particle conversion (number).
        nprci(i) = n0i(i)/(lami(i)*ac_time(i))*exp(-d_rat)

        m_ip = rhoi * pio6 / (lami(i)*lami(i)*lami(i))

!       m_ip = (rhoi*pi/6._r8) / lami(i)**3

        ! Rate of mass conversion.
        ! Note that this is:
        ! m n (d^3 + 3 d^2 + 6 d + 6)
        prci(i) = m_ip * nprci(i) * (((d_rat + three)*d_rat + six)*d_rat + six)

     else
        prci(i)  = zero
        nprci(i) = zero
     end if
  enddo
end subroutine ice_autoconversion
!===================================
! Anning Cheng 10/5/2017 added GMAO ice autoconversion
!>\ingroup micro_mg_utils_mod
!! GMAO ice autoconversion
subroutine gmao_ice_autoconversion(t, qiic, niic, lami, n0i,        &
                                   dcs, ac_time, prci, nprci, mgncol)

      integer, intent(in) :: mgncol
      real(r8), dimension(mgncol), intent(in) :: t
      real(r8), dimension(mgncol), intent(in) :: qiic
      real(r8), dimension(mgncol), intent(in) :: niic
      real(r8), dimension(mgncol), intent(in) :: lami
      real(r8), dimension(mgncol), intent(in) :: n0i
      real(r8), dimension(mgncol), intent(in) :: ac_time 
      real(r8), intent(in) :: dcs

      real(r8), dimension(mgncol), intent(out) :: prci
      real(r8), dimension(mgncol), intent(out) :: nprci


      real(r8) :: m_ip, tx1, tx2
      integer  :: i
      do i=1,mgncol
        if (t(i) <= tmelt .and. qiic(i) >= qsmall) then
           m_ip = max(min(0.008_r8*(lami(i)*0.01)**0.87_r8,  &
                          10.0_r8), 0.1_r8)
           tx1 = lami(i)*dcs
           tx2 = one / ac_time(i)
           nprci(i) = niic(i)*tx2 * (one - gamma_incomp(m_ip, tx1))
           prci(i)  = qiic(i)*tx2 * (one - gamma_incomp(m_ip+three, tx1))
         else
           prci(i)  = zero
           nprci(i) = zero
         end if
       enddo
end subroutine gmao_ice_autoconversion
!===================================
! immersion freezing (Bigg, 1953)
!===================================
!>\ingroup micro_mg_utils_mod
!! immersion freezing (Bigg, 1953)
subroutine immersion_freezing(microp_uniform, t, pgam, lamc, &
                              qcic, ncic, relvar, mnuccc, nnuccc, mgncol)

  integer, intent(in) :: mgncol
  logical, intent(in) :: microp_uniform

  ! Temperature
  real(r8), dimension(mgncol), intent(in) :: t

  ! Cloud droplet size distribution parameters
  real(r8), dimension(mgncol), intent(in) :: pgam
  real(r8), dimension(mgncol), intent(in) :: lamc

  ! MMR and number concentration of in-cloud liquid water
  real(r8), dimension(mgncol), intent(in) :: qcic
  real(r8), dimension(mgncol), intent(in) :: ncic

  ! Relative variance of cloud water
  real(r8), dimension(mgncol), intent(in) :: relvar

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: mnuccc ! MMR
  real(r8), dimension(mgncol), intent(out) :: nnuccc ! Number

  ! Coefficients that will be omitted for sub-columns
  real(r8), dimension(mgncol) :: dum
  real(r8)                    :: tx1
  integer :: i

  if (.not. microp_uniform) then
     dum(:) = var_coef(relvar, 2)
  else
     dum(:) = one
  end if
  do i=1,mgncol

     if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then

        tx1 = one / (lamc(i) * lamc(i) * lamc(i))
        nnuccc(i) = pio6*ncic(i)*rising_factorial(pgam(i)+one, 3) * &
                     bimm*(exp(aimm*(tmelt - t(i)))-one) * tx1

        mnuccc(i) = dum(i) * nnuccc(i) * pio6 * rhow * &
                    rising_factorial(pgam(i)+four, 3) * tx1

     else
        mnuccc(i) = zero
        nnuccc(i) = zero
     end if ! qcic > qsmall and t < 4 deg C
  enddo

end subroutine immersion_freezing

!>\ingroup micro_mg_utils_mod
!! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
!! dust size and number in multiple bins are read in from companion routine
subroutine contact_freezing (microp_uniform, t, p, rndst, nacon, &
     pgam, lamc, qcic, ncic, relvar, mnucct, nnucct, mgncol, mdust)

  logical, intent(in) :: microp_uniform

  integer, intent(in) :: mgncol
  integer, intent(in) :: mdust

  real(r8), dimension(mgncol), intent(in) :: t            ! Temperature
  real(r8), dimension(mgncol), intent(in) :: p            ! Pressure
  real(r8), dimension(mgncol, mdust), intent(in) :: rndst ! Radius (for multiple dust bins)
  real(r8), dimension(mgncol, mdust), intent(in) :: nacon ! Number (for multiple dust bins)

  ! Size distribution parameters for cloud droplets
  real(r8), dimension(mgncol), intent(in) :: pgam
  real(r8), dimension(mgncol), intent(in) :: lamc

  ! MMR and number concentration of in-cloud liquid water
  real(r8), dimension(mgncol), intent(in) :: qcic
  real(r8), dimension(mgncol), intent(in) :: ncic

  ! Relative cloud water variance
  real(r8), dimension(mgncol), intent(in) :: relvar

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: mnucct ! MMR
  real(r8), dimension(mgncol), intent(out) :: nnucct ! Number

  real(r8) :: tcnt                  ! scaled relative temperature
  real(r8) :: viscosity             ! temperature-specific viscosity (kg/m/s)
  real(r8) :: mfp                   ! temperature-specific mean free path (m)

  ! Dimension these according to number of dust bins, inferred from rndst size
  real(r8) :: nslip(size(rndst,2))  ! slip correction factors
  real(r8) :: ndfaer(size(rndst,2)) ! aerosol diffusivities (m^2/sec)

  ! Coefficients not used for subcolumns
  real(r8) :: dum, dum1, tx1

  ! Common factor between mass and number.
  real(r8) :: contact_factor

  integer  :: i

  do i = 1,mgncol

     if (qcic(i) >= qsmall .and. t(i) < 269.15_r8) then

        if (.not. microp_uniform) then
           dum  = var_coef(relvar(i), four/three)
           dum1 = var_coef(relvar(i), oneo3)
        else
           dum  = one
           dum1 = one
        endif

        tcnt=(270.16_r8-t(i))**1.3_r8
        viscosity = 1.8e-5_r8*(t(i)/298.0_r8)**0.85_r8    ! Viscosity (kg/m/s)
        mfp = two*viscosity/ &                         ! Mean free path (m)
                     (p(i)*sqrt( 8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i)) ))

        ! Note that these two are vectors.
        nslip = one+(mfp/rndst(i,:))*(1.257_r8+(0.4_r8*exp(-(1.1_r8*rndst(i,:)/mfp))))! Slip correction factor

        ndfaer = 1.381e-23_r8*t(i)*nslip/(6._r8*pi*viscosity*rndst(i,:))  ! aerosol diffusivity (m2/s)

        tx1 = one / lamc(i)
        contact_factor = dot_product(ndfaer,nacon(i,:)*tcnt) * pi * &
                         ncic(i) * (pgam(i) + one) * tx1

        mnucct(i) = dum * contact_factor * &
                    pio3*rhow*rising_factorial(pgam(i)+two, 3) * tx1 * tx1 *tx1

        nnucct(i) =  (dum1+dum1) * contact_factor

     else

        mnucct(i) = zero
        nnucct(i) = zero

     end if ! qcic > qsmall and t < 4 deg C
  end do

end subroutine contact_freezing

!>\ingroup micro_mg_utils_mod
!! snow self-aggregation from passarelli, 1978, used by reisner, 1998
!===================================================================
! this is hard-wired for bs = 0.4 for now
! ignore self-collection of cloud ice
subroutine snow_self_aggregation(t, rho, asn, rhosn, qsic, nsic, nsagg, mgncol)

  integer,                          intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t     ! Temperature
  real(r8), dimension(mgncol), intent(in) :: rho   ! Density
  real(r8), dimension(mgncol), intent(in) :: asn   ! fall speed parameter for snow
  real(r8),                    intent(in) :: rhosn ! density of snow

  ! In-cloud snow
  real(r8), dimension(mgncol), intent(in) :: qsic ! MMR
  real(r8), dimension(mgncol), intent(in) :: nsic ! Number

  ! Output number tendency
  real(r8), dimension(mgncol), intent(out) :: nsagg

  integer :: i

  do i=1,mgncol
     if (qsic(i) >= qsmall .and. t(i) <= tmelt) then
        nsagg(i) = -1108._r8*eii/(four*720._r8*rhosn)*asn(i)*qsic(i)*nsic(i)*rho(i)*&
             ((qsic(i)/nsic(i))*(one/(rhosn*pi)))**((bs-one)*oneo3)
     else
        nsagg(i) = zero
     end if
  enddo
end subroutine snow_self_aggregation

!>\ingroup micro_mg_utils_mod
!! accretion of cloud droplets onto snow/graupel
!===================================================================
! here use continuous collection equation with
! simple gravitational collection kernel
! ignore collisions between droplets/cloud ice
! since minimum size ice particle for accretion is 50 - 150 micron
subroutine accrete_cloud_water_snow(t, rho, asn, uns, mu, qcic, ncic, qsic, &
     pgam, lamc, lams, n0s, psacws, npsacws, mgncol)

  integer, intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t   ! Temperature
  real(r8), dimension(mgncol), intent(in) :: rho ! Density
  real(r8), dimension(mgncol), intent(in) :: asn ! Fallspeed parameter (snow)
  real(r8), dimension(mgncol), intent(in) :: uns ! Current fallspeed   (snow)
  real(r8), dimension(mgncol), intent(in) :: mu  ! Viscosity

  ! In-cloud liquid water
  real(r8), dimension(mgncol), intent(in) :: qcic ! MMR
  real(r8), dimension(mgncol), intent(in) :: ncic ! Number

  ! In-cloud snow
  real(r8), dimension(mgncol), intent(in) :: qsic ! MMR

  ! Cloud droplet size parameters
  real(r8), dimension(mgncol), intent(in) :: pgam
  real(r8), dimension(mgncol), intent(in) :: lamc

  ! Snow size parameters
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: psacws  ! Mass mixing ratio
  real(r8), dimension(mgncol), intent(out) :: npsacws ! Number concentration

  real(r8) :: dc0 ! Provisional mean droplet size
  real(r8) :: dum
  real(r8) :: eci ! collection efficiency for riming of snow by droplets

  ! Fraction of cloud droplets accreted per second
  real(r8) :: accrete_rate
  integer :: i

  ! ignore collision of snow with droplets above freezing

  do i=1,mgncol
     if (qsic(i) >= qsmall .and. t(i) <= tmelt .and. qcic(i) >= qsmall) then

        ! put in size dependent collection efficiency
        ! mean diameter of snow is area-weighted, since
        ! accretion is function of crystal geometric area
        ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)

        dc0 = (pgam(i)+one)/lamc(i)
        dum = dc0*dc0*uns(i)*rhow*lams(i)/(9._r8*mu(i))
        eci = dum*dum / ((dum+0.4_r8)*(dum+0.4_r8))

        eci = max(eci,zero)
        eci = min(eci,one)

        ! no impact of sub-grid distribution of qc since psacws
        ! is linear in qc
        accrete_rate = (pi/four)*asn(i)*rho(i)*n0s(i)*eci*gamma_bs_plus3 / lams(i)**(bs+three)
        psacws(i)    = accrete_rate*qcic(i)
        npsacws(i)   = accrete_rate*ncic(i)
     else
        psacws(i)  = zero
        npsacws(i) = zero
     end if
  enddo
end subroutine accrete_cloud_water_snow

!>\ingroup micro_mg_utils_mod
!! add secondary ice production due to accretion of droplets by snow
!===================================================================
! (Hallet-Mossop process) (from Cotton et al., 1986)
subroutine secondary_ice_production(t, psacws, msacwi, nsacwi, mgncol)

  integer, intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t ! Temperature

  ! Accretion of cloud water to snow tendencies
  real(r8), dimension(mgncol), intent(inout) :: psacws ! MMR

  ! Output (ice) tendencies
  real(r8), dimension(mgncol), intent(out) :: msacwi ! MMR
  real(r8), dimension(mgncol), intent(out) :: nsacwi ! Number
  integer :: i

  do i=1,mgncol
     if((t(i) < 270.16_r8) .and. (t(i) >= 268.16_r8)) then
        nsacwi(i) = 3.5e8_r8*(270.16_r8-t(i))/two*psacws(i)
     else if((t(i) < 268.16_r8) .and. (t(i) >= 265.16_r8)) then
        nsacwi(i) = 3.5e8_r8*(t(i)-265.16_r8)*oneo3*psacws(i)
     else
        nsacwi(i) = zero
     endif
  enddo

  do i=1,mgncol
     msacwi(i) = min(nsacwi(i)*mi0, psacws(i))
     psacws(i) = psacws(i) - msacwi(i)
  enddo
end subroutine secondary_ice_production

!>\ingroup micro_mg_utils_mod
!! accretion of rain water by snow
!===================================================================
! formula from ikawa and saito, 1991, used by reisner et al., 1998
subroutine accrete_rain_snow(t, rho, umr, ums, unr, uns, qric, qsic, &
     lamr, n0r, lams, n0s, pracs, npracs, mgncol)

  integer,                          intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t   ! Temperature
  real(r8), dimension(mgncol), intent(in) :: rho ! Density

  ! Fallspeeds
  ! mass-weighted
  real(r8), dimension(mgncol), intent(in) :: umr ! rain
  real(r8), dimension(mgncol), intent(in) :: ums ! snow
  ! number-weighted
  real(r8), dimension(mgncol), intent(in) :: unr ! rain
  real(r8), dimension(mgncol), intent(in) :: uns ! snow

  ! In cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qric ! rain
  real(r8), dimension(mgncol), intent(in) :: qsic ! snow

  ! Size distribution parameters
  ! rain
  real(r8), dimension(mgncol), intent(in) :: lamr
  real(r8), dimension(mgncol), intent(in) :: n0r
  ! snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: pracs  ! MMR
  real(r8), dimension(mgncol), intent(out) :: npracs ! Number

  ! Collection efficiency for accretion of rain by snow
  real(r8), parameter :: ecr = one

  ! Ratio of average snow diameter to average rain diameter.
  real(r8) :: d_rat
  ! Common factor between mass and number expressions
  real(r8) :: common_factor
  real(r8) :: tx1, tx2
  integer :: i

  do i=1,mgncol
     if (qric(i) >= icsmall .and. qsic(i) >= icsmall .and. t(i) <= tmelt) then

        tx2 = lamr(i)*lamr(i)*lamr(i)

        common_factor = pi*ecr*rho(i)*n0r(i)*n0s(i) / (tx2 * lams(i))

        d_rat = lamr(i)/lams(i)

        tx1       = 1.2_r8*umr(i)-0.95_r8*ums(i)
        pracs(i)  = common_factor*pi*rhow*                  &
                    sqrt(tx1*tx1 + 0.08_r8*ums(i)*umr(i)) * &
                   ((half*d_rat + two)*d_rat + five) / tx2

        tx1       = unr(i)-uns(i)
        npracs(i) = common_factor*half *                          &
                    sqrt(1.7_r8*tx1*tx1 + 0.3_r8*unr(i)*uns(i)) * &
                   ((d_rat + one)*d_rat + one)

     else
        pracs(i)  = zero
        npracs(i) = zero
     end if
  enddo
end subroutine accrete_rain_snow

!>\ingroup micro_mg_utils_mod
!! heterogeneous freezing of rain drops
!===================================================================
! follows from Bigg (1953)
subroutine heterogeneous_rain_freezing(t, qric, nric, lamr, mnuccr, nnuccr, mgncol)

  integer,                          intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t    ! Temperature

  ! In-cloud rain
  real(r8), dimension(mgncol), intent(in) :: qric ! MMR
  real(r8), dimension(mgncol), intent(in) :: nric ! Number
  real(r8), dimension(mgncol), intent(in) :: lamr ! size parameter

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: mnuccr ! MMR
  real(r8), dimension(mgncol), intent(out) :: nnuccr ! Number
  real(r8) :: tx1
  integer  :: i

  do i=1,mgncol

     if (t(i) < 269.15_r8 .and. qric(i) >= qsmall) then
        tx1 = pi / (lamr(i)*lamr(i)*lamr(i))
        nnuccr(i) = nric(i)*bimm* (exp(aimm*(tmelt - t(i)))-one) * tx1

        mnuccr(i) = nnuccr(i) * 20._r8*rhow * tx1

     else
        mnuccr(i) = zero
        nnuccr(i) = zero
     end if
  enddo
end subroutine heterogeneous_rain_freezing

!>\ingroup micro_mg_utils_mod
!! accretion of cloud liquid water by rain
!! formula from Khrouditnov and Kogan (2000)
! gravitational collection kernel, droplet fall speed neglected
subroutine accrete_cloud_water_rain(microp_uniform, qric, qcic, &
     ncic, relvar, accre_enhan, pra, npra, mgncol)

  logical, intent(in) :: microp_uniform
  integer, intent(in) :: mgncol
  ! In-cloud rain
  real(r8), dimension(mgncol), intent(in) :: qric ! MMR

  ! Cloud droplets
  real(r8), dimension(mgncol), intent(in) :: qcic ! MMR
  real(r8), dimension(mgncol), intent(in) :: ncic ! Number

  ! SGS variability
  real(r8), dimension(mgncol), intent(in) :: relvar
  real(r8), dimension(mgncol), intent(in) :: accre_enhan

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: pra  ! MMR
  real(r8), dimension(mgncol), intent(out) :: npra ! Number

  ! Coefficient that varies for subcolumns
  real(r8), dimension(mgncol) :: pra_coef

  integer :: i

  if (.not. microp_uniform) then
    pra_coef(:) = accre_enhan * var_coef(relvar(:), 1.15_r8)
  else
    pra_coef(:) = one
  end if

  do i=1,mgncol

    if (qric(i) >= qsmall .and. qcic(i) >= qsmall) then

      ! include sub-grid distribution of cloud water
      pra(i) = pra_coef(i) * 67._r8*(qcic(i)*qric(i))**1.15_r8

      npra(i) = pra(i)*ncic(i)/qcic(i)

    else
      pra(i)  = zero
      npra(i) = zero
    end if
  end do
end subroutine accrete_cloud_water_rain

!>\ingroup micro_mg_utils_mod
!! Self-collection of rain drops
!! from Beheng(1994)
subroutine self_collection_rain(rho, qric, nric, nragg, mgncol)

  integer,                          intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: rho  ! Air density

  ! Rain
  real(r8), dimension(mgncol), intent(in) :: qric ! MMR
  real(r8), dimension(mgncol), intent(in) :: nric ! Number

  ! Output number tendency
  real(r8), dimension(mgncol), intent(out) :: nragg

  integer :: i

  do i=1,mgncol
     if (qric(i) >= qsmall) then
        nragg(i) = -8._r8*nric(i)*qric(i)*rho(i)
     else
        nragg(i) = zero
     end if
  enddo
end subroutine self_collection_rain

!>\ingroup micro_mg_utils_mod
!! Accretion of cloud ice by snow
!===================================================================
! For this calculation, it is assumed that the Vs >> Vi
! and Ds >> Di for continuous collection
subroutine accrete_cloud_ice_snow(t, rho, asn, qiic, niic, qsic, &
     lams, n0s, prai, nprai, mgncol)

  integer,                          intent(in) :: mgncol
  real(r8), dimension(mgncol), intent(in) :: t    ! Temperature
  real(r8), dimension(mgncol), intent(in) :: rho   ! Density

  real(r8), dimension(mgncol), intent(in) :: asn  ! Snow fallspeed parameter

  ! Cloud ice
  real(r8), dimension(mgncol), intent(in) :: qiic ! MMR
  real(r8), dimension(mgncol), intent(in) :: niic ! Number

  real(r8), dimension(mgncol), intent(in) :: qsic ! Snow MMR

  ! Snow size parameters
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: prai ! MMR
  real(r8), dimension(mgncol), intent(out) :: nprai ! Number

  ! Fraction of cloud ice particles accreted per second
  real(r8) :: accrete_rate

  integer :: i

  do i=1,mgncol
     if (qsic(i) >= qsmall .and. qiic(i) >= qsmall .and. t(i) <= tmelt) then

        accrete_rate = (pi/four) * eii * asn(i) * rho(i) * n0s(i) * gamma_bs_plus3  &
                     / lams(i)**(bs+three)

        prai(i)  = accrete_rate * qiic(i)
        nprai(i) = accrete_rate * niic(i)

     else
        prai(i)  = zero
        nprai(i) = zero
     end if
  enddo
end subroutine accrete_cloud_ice_snow

!>\ingroup micro_mg_utils_mod
!! calculate evaporation/sublimation of rain and snow
!===================================================================
! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
! in-cloud condensation/deposition of rain and snow is neglected
! except for transfer of cloud water to snow through bergeron process
subroutine evaporate_sublimate_precip(t, rho, dv, mu, sc, q, qvl, qvi, &
     lcldm, precip_frac, arn, asn, qcic, qiic, qric, qsic, lamr, n0r, lams, n0s, &
     pre, prds, am_evp_st, mgncol)

  integer,  intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t    ! temperature
  real(r8), dimension(mgncol), intent(in) :: rho  ! air density
  real(r8), dimension(mgncol), intent(in) :: dv   ! water vapor diffusivity
  real(r8), dimension(mgncol), intent(in) :: mu   ! viscosity
  real(r8), dimension(mgncol), intent(in) :: sc   ! schmidt number
  real(r8), dimension(mgncol), intent(in) :: q    ! humidity
  real(r8), dimension(mgncol), intent(in) :: qvl  ! saturation humidity (water)
  real(r8), dimension(mgncol), intent(in) :: qvi  ! saturation humidity (ice)
  real(r8), dimension(mgncol), intent(in) :: lcldm  ! liquid cloud fraction
  real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)

  ! fallspeed parameters
  real(r8), dimension(mgncol), intent(in) :: arn  ! rain
  real(r8), dimension(mgncol), intent(in) :: asn  ! snow

  ! In-cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid
  real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice
  real(r8), dimension(mgncol), intent(in) :: qric ! rain
  real(r8), dimension(mgncol), intent(in) :: qsic ! snow

  ! Size parameters
  ! rain
  real(r8), dimension(mgncol), intent(in) :: lamr
  real(r8), dimension(mgncol), intent(in) :: n0r
  ! snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: pre
  real(r8), dimension(mgncol), intent(out) :: prds
  real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates.

  real(r8) :: qclr   ! water vapor mixing ratio in clear air
  real(r8) :: ab     ! correction to account for latent heat
  real(r8) :: eps    ! 1/ sat relaxation timescale
  real(r8) :: tx1, tx2, tx3

  real(r8), dimension(mgncol) :: dum

  integer :: i

  am_evp_st = zero
  ! set temporary cloud fraction to zero if cloud water + ice is very small
  ! this will ensure that evaporation/sublimation of precip occurs over
  ! entire grid cell, since min cloud fraction is specified otherwise
  do i=1,mgncol
     if (qcic(i)+qiic(i) < 1.e-6_r8) then
        dum(i) = zero
     else
        dum(i) = lcldm(i)
     end if
  enddo
  do i=1,mgncol
  ! only calculate if there is some precip fraction > cloud fraction

     if (precip_frac(i) > dum(i)) then

        if (qric(i) >= qsmall .or. qsic(i) >= qsmall) then
           am_evp_st(i) = precip_frac(i) - dum(i)

           ! calculate q for out-of-cloud region
           qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i))
        end if

        ! evaporation of rain
        if (qric(i) >= qsmall) then

           ab  = calc_ab(t(i), qvl(i), xxlv)
           eps = two*pi*n0r(i)*rho(i)*Dv(i) *             &
                (f1r/(lamr(i)*lamr(i)) +                  &
                 f2r*sqrt(arn(i)*rho(i)/mu(i)) *          &
                 sc(i)**oneo3*gamma_half_br_plus5         &
                 / (lamr(i)**((five+br)*half)))

           pre(i) = eps*(qclr-qvl(i)) / ab

           ! only evaporate in out-of-cloud region
           ! and distribute across precip_frac
           pre(i) = min(pre(i)*am_evp_st(i), zero)
           pre(i) = pre(i) / precip_frac(i)
        else
           pre(i) = zero
        end if

        ! sublimation of snow
        if (qsic(i) >= qsmall) then
           ab = calc_ab(t(i), qvi(i), xxls)
           eps = two*pi*n0s(i)*rho(i)*Dv(i) *             &
                ( f1s/(lams(i)*lams(i))                   &
                + f2s*sqrt(asn(i)*rho(i)/mu(i)) *         &
                  sc(i)**oneo3*gamma_half_bs_plus5        &
                / (lams(i)**((five+bs)*half)))
           prds(i) = eps*(qclr-qvi(i)) / ab

           ! only sublimate in out-of-cloud region and distribute over precip_frac
           prds(i) = min(prds(i)*am_evp_st(i), zero)
           prds(i) = prds(i) / precip_frac(i)
        else
           prds(i) = zero
        end if

     else
        prds(i) = zero
        pre(i)  = zero
     end if
  enddo

end subroutine evaporate_sublimate_precip

!>\ingroup micro_mg_utils_mod
!! evaporation/sublimation of rain, snow and graupel
!===================================================================
! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
! in-cloud condensation/deposition of rain and snow is neglected
! except for transfer of cloud water to snow through bergeron process
subroutine evaporate_sublimate_precip_graupel(t, rho, dv, mu, sc, q, qvl, qvi, &
     lcldm, precip_frac, arn, asn, agn, bg, qcic, qiic, qric, qsic, qgic, lamr, n0r, lams, n0s, lamg, n0g, &
     pre, prds, prdg, am_evp_st, mgncol)

  integer,  intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t    ! temperature
  real(r8), dimension(mgncol), intent(in) :: rho  ! air density
  real(r8), dimension(mgncol), intent(in) :: dv   ! water vapor diffusivity
  real(r8), dimension(mgncol), intent(in) :: mu   ! viscosity
  real(r8), dimension(mgncol), intent(in) :: sc   ! schmidt number
  real(r8), dimension(mgncol), intent(in) :: q    ! humidity
  real(r8), dimension(mgncol), intent(in) :: qvl  ! saturation humidity (water)
  real(r8), dimension(mgncol), intent(in) :: qvi  ! saturation humidity (ice)
  real(r8), dimension(mgncol), intent(in) :: lcldm  ! liquid cloud fraction
  real(r8), dimension(mgncol), intent(in) :: precip_frac ! precipitation fraction (maximum overlap)

  ! fallspeed parameters
  real(r8), dimension(mgncol), intent(in) :: arn  ! rain
  real(r8), dimension(mgncol), intent(in) :: asn  ! snow
!++ag
  real(r8), dimension(mgncol), intent(in) :: agn  ! graupel
  real(r8),                    intent(in) :: bg
!--ag

  ! In-cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid
  real(r8), dimension(mgncol), intent(in) :: qiic ! cloud ice
  real(r8), dimension(mgncol), intent(in) :: qric ! rain
  real(r8), dimension(mgncol), intent(in) :: qsic ! snow
  real(r8), dimension(mgncol), intent(in) :: qgic ! graupel

  ! Size parameters
  ! rain
  real(r8), dimension(mgncol), intent(in) :: lamr
  real(r8), dimension(mgncol), intent(in) :: n0r
  ! snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s
!++ag
  ! graupel
  real(r8), dimension(mgncol), intent(in) :: lamg
  real(r8), dimension(mgncol), intent(in) :: n0g
!--ag

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: pre
  real(r8), dimension(mgncol), intent(out) :: prds
!++ag
  real(r8), dimension(mgncol), intent(out) :: prdg
!--ag
  real(r8), dimension(mgncol), intent(out) :: am_evp_st ! Fractional area where rain evaporates.

  real(r8) :: qclr   ! water vapor mixing ratio in clear air
  real(r8) :: ab     ! correction to account for latent heat
  real(r8) :: eps    ! 1/ sat relaxation timescale

  real(r8), dimension(mgncol) :: dum

  integer :: i

  ! set temporary cloud fraction to zero if cloud water + ice is very small
  ! this will ensure that evaporation/sublimation of precip occurs over
  ! entire grid cell, since min cloud fraction is specified otherwise
  do i=1,mgncol
     if (qcic(i)+qiic(i) < 1.e-6_r8) then
        dum(i) = zero
     else
        dum(i) = lcldm(i)
     end if
  enddo
  do i=1,mgncol
  ! only calculate if there is some precip fraction > cloud fraction

     if (precip_frac(i) > dum(i)) then

        if (qric(i) >= qsmall .or. qsic(i) >= qsmall .or. qgic(i) >= qsmall) then
           am_evp_st(i) = precip_frac(i) - dum(i)

           ! calculate q for out-of-cloud region
           qclr = (q(i)-dum(i)*qvl(i)) / (one-dum(i))
        end if

        ! evaporation of rain
        if (qric(i) >= qsmall) then

           ab  = calc_ab(t(i), qvl(i), xxlv)
           eps = twopi*n0r(i)*rho(i)*Dv(i)*        &
                ( f1r/(lamr(i)*lamr(i))            &
                + f2r*sqrt(arn(i)*rho(i)/mu(i))    &
                * sc(i)**oneo3*gamma_half_br_plus5 &
                / (lamr(i)**((five+br)*half)))

           pre(i) = eps*(qclr-qvl(i))/ab

           ! only evaporate in out-of-cloud region
           ! and distribute across precip_frac
           pre(i) = min(pre(i)*am_evp_st(i), zero)
           pre(i) = pre(i)/precip_frac(i)
        else
           pre(i) = zero
        end if

        ! sublimation of snow
        if (qsic(i) >= qsmall) then
           ab  = calc_ab(t(i), qvi(i), xxls)
           eps = twopi*n0s(i)*rho(i)*Dv(i)*        &
                ( f1s/(lams(i)*lams(i))            &
                + f2s*sqrt(asn(i)*rho(i)/mu(i))    &
                * sc(i)**oneo3*gamma_half_bs_plus5 &
                / (lams(i)**((five+bs)*half)))
           prds(i) = eps*(qclr-qvi(i))/ab

          ! only sublimate in out-of-cloud region and distribute over precip_frac
           prds(i) = min(prds(i)*am_evp_st(i), zero)
           prds(i) = prds(i)/precip_frac(i)
        else
           prds(i) = zero
        end if

!++AG ADD GRAUPEL, do Same with prdg.

        if (qgic(i) >= qsmall) then
           ab = calc_ab(t(i), qvi(i), xxls)

           eps = twopi*n0g(i)*rho(i)*Dv(i)*          &
                ( f1s/(lamg(i)*lamg(i))              &
                + f2s*sqrt(agn(i)*rho(i)/mu(i))      &
                * sc(i)**oneo3*gamma((five+bg)*half) &
                / (lamg(i)**((five+bs)*half)))
!               / (lamg(i)**((five+bg)*half)))              ! changing bs to bg - Moorthi
           prdg(i) = eps*(qclr-qvi(i))/ab

           ! only sublimate in out-of-cloud region and distribute over precip_frac
           prdg(i) = min(prdg(i)*am_evp_st(i), zero)
           prdg(i) = prdg(i)/precip_frac(i)
        else
           prdg(i) = zero
        end if

     else
        prds(i) = zero
        pre(i)  =  zero
!++ag
        prdg(i) = zero
!--ag
     end if
  enddo

end subroutine evaporate_sublimate_precip_graupel

!>\ingroup micro_mg_utils_mod
!! bergeron process - evaporation of droplets and deposition onto snow
subroutine bergeron_process_snow(t, rho, dv, mu, sc, qvl, qvi, asn, &
                                 qcic, qsic, lams, n0s, bergs, mgncol)

  integer, intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t    ! temperature
  real(r8), dimension(mgncol), intent(in) :: rho  ! air density
  real(r8), dimension(mgncol), intent(in) :: dv   ! water vapor diffusivity
  real(r8), dimension(mgncol), intent(in) :: mu   ! viscosity
  real(r8), dimension(mgncol), intent(in) :: sc   ! schmidt number
  real(r8), dimension(mgncol), intent(in) :: qvl  ! saturation humidity (water)
  real(r8), dimension(mgncol), intent(in) :: qvi  ! saturation humidity (ice)

  ! fallspeed parameter for snow
  real(r8), dimension(mgncol), intent(in) :: asn

  ! In-cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio

  ! Size parameters for snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: bergs

  real(r8) :: ab     ! correction to account for latent heat
  real(r8) :: eps    ! 1/ sat relaxation timescale

  integer :: i

  do i=1,mgncol
     if (qsic(i) >= qsmall.and. qcic(i) >= qsmall .and. t(i) < tmelt) then
        ab = calc_ab(t(i), qvi(i), xxls)
        eps = two*pi*n0s(i)*rho(i)*Dv(i) *            &
             (f1s/(lams(i)*lams(i)) +                 &
             f2s*sqrt(asn(i)*rho(i)/mu(i)) *          &
             sc(i)**oneo3*gamma_half_bs_plus5 /       &
             (lams(i)**((five+bs)*half)))
        bergs(i) = eps*(qvl(i)-qvi(i)) / ab
     else
        bergs(i) = zero
     end if
  enddo
end subroutine bergeron_process_snow

!========================================================================
!>\ingroup micro_mg_utils_mod
!! Collection of snow by rain to form graupel
subroutine graupel_collecting_snow(qsic,qric,umr,ums,rho,lamr,n0r,lams,n0s, &
                                   psacr, mgncol)

  integer, intent(in) :: mgncol

  ! In-cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qsic ! snow
  real(r8), dimension(mgncol), intent(in) :: qric ! rain

  ! mass-weighted fall speeds
  real(r8), dimension(mgncol), intent(in) :: umr ! rain
  real(r8), dimension(mgncol), intent(in) :: ums ! snow

  real(r8), dimension(mgncol), intent(in) :: rho  ! air density


  ! Size parameters for rain
  real(r8), dimension(mgncol), intent(in) :: lamr
  real(r8), dimension(mgncol), intent(in) :: n0r

  ! Size parameters for snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  real(r8), dimension(mgncol), intent(out) :: psacr ! conversion due to coll of snow by rain

  real(r8) :: cons31, tx1, tx2, tx3, tx4, tx5
  integer :: i

  cons31 = pi*pi*ecr*rhosn

  do i=1,mgncol

     if (qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then
        tx1 = 1.2_r8*umr(i) - 0.95_r8*ums(i)
        tx1 = sqrt(tx1*tx1+0.08_r8*ums(i)*umr(i))
        tx2 = one / lams(i)
        tx3 = one / lamr(i)
        tx4 = tx2 * tx2
        tx5 = tx4 * tx4 * tx3

        psacr(i) = cons31 * tx1 * rho(i) * n0r(i) * n0s(i) * tx5    &
                                * (5.0_r8*tx4+tx3*(tx2+tx2+0.5_r8*tx3))

!       psacr(i) = cons31*(((1.2_r8*umr(i)-0.95_r8*ums(i))**2+      &
!            0.08_r8*ums(i)*umr(i))**0.5_r8*rho(i)*                 &
!            n0r(i)*n0s(i)/lams(i)**3*                              &
!            (5._r8/(lams(i)**3*lamr(i))+                           &
!            2._r8/(lams(i)**2*lamr(i)**2)+                         &
!            0.5_r8/(lams(i)*lamr(i)**3)))
     else
        psacr(i) = zero
     end if

  end do

end subroutine graupel_collecting_snow

!========================================================================
!>\ingroup micro_mg_utils_mod
!! Collection of cloud water by graupel
subroutine graupel_collecting_cld_water(qgic,qcic,ncic,rho,n0g,lamg,bg,agn, &
                                        psacwg, npsacwg, mgncol)

  integer, intent(in) :: mgncol

  ! In-cloud MMRs
  real(r8), dimension(mgncol), intent(in) :: qgic ! graupel
  real(r8), dimension(mgncol), intent(in) :: qcic ! cloud water

  real(r8), dimension(mgncol), intent(in) :: ncic ! cloud water number conc

  real(r8), dimension(mgncol), intent(in) :: rho  ! air density

  ! Size parameters for graupel
  real(r8), dimension(mgncol), intent(in) :: lamg
  real(r8), dimension(mgncol), intent(in) :: n0g

  ! fallspeed parameters for graupel
  ! Set AGN and BG  as input (in micro_mg3_0.F90)  (select hail or graupel as appropriate)
  real(r8),                    intent(in) :: bg
  real(r8), dimension(mgncol), intent(in) :: agn

  ! Output tendencies
  real(r8), dimension(mgncol), intent(out) :: psacwg
  real(r8), dimension(mgncol), intent(out) :: npsacwg

  real(r8) :: cons, tx1
  integer :: i

  cons = gamma(bg+three) * pi/four * ecid

  do i=1,mgncol

        if (qgic(i) >= 1.e-8_r8 .and. qcic(i) >= qsmall) then

           tx1 = cons*agn(i)*rho(i)*n0g(i) / lamg(i)**(bg+three)

           psacwg(i)  = tx1 * qcic(i)
           npsacwg(i) = tx1 * ncic(i)
        else
           psacwg(i)  = zero
           npsacwg(i) = zero
        end if
     enddo
end subroutine graupel_collecting_cld_water

!========================================================================
!>\ingroup micro_mg_utils_mod
!! Conversion of rimed cloud water onto snow to graupel/hail
subroutine graupel_riming_liquid_snow(psacws,qsic,qcic,nsic,rho,rhosn,rhog,asn,lams,n0s,dtime, &
                                      pgsacw,nscng,mgncol)

  integer, intent(in) :: mgncol

  ! Accretion of cloud water to snow tendency (modified)
  real(r8), dimension(mgncol), intent(inout) :: psacws

  real(r8), dimension(mgncol), intent(in) :: qsic ! snow mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qcic ! cloud liquid mixing ratio
  real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration

  real(r8), dimension(mgncol), intent(in) :: rho   ! air density
  real(r8),                    intent(in) :: rhosn ! snow density
  real(r8),                    intent(in) :: rhog ! graupel density

  real(r8), dimension(mgncol), intent(in) :: asn   ! fall speed parameter for snow

  ! Size parameters for snow
  real(r8), dimension(mgncol), intent(in) :: lams
  real(r8), dimension(mgncol), intent(in) :: n0s

  real(r8),                    intent(in) :: dtime

  !Output process rates
  real(r8), dimension(mgncol), intent(out) :: pgsacw  ! dQ graupel due to collection droplets by snow
  real(r8), dimension(mgncol), intent(out) :: nscng   ! dN graupel due to collection droplets by snow

  real(r8) :: cons
  real(r8) :: rhosu
  real(r8) :: dum
  integer :: i

!........................................................................
!Input: PSACWS,qs,qc,n0s,rho,lams,rhos,rhog
!Output:PSACWS,PGSACW,NSCNG

  rhosu = 85000._r8/(ra * tmelt)    ! typical air density at 850 mb

  do i=1,mgncol

       cons=4._r8 *2._r8 *3._r8 *rhosu*pi*ecid*ecid*gamma_2bs_plus2/(8._r8*(rhog-rhosn))

     if (psacws(i).gt.0._r8 .and. qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then
! Only allow conversion if qni > 0.1 and qc > 0.5 g/kg following Rutledge and Hobbs (1984)
        !if (qsic(i).GE.0.1e-3_r8 .AND. qcic(i).GE.0.5E-3_r8) then

! portion of riming converted to graupel (Reisner et al. 1998, originally IS1991)
! dtime here is correct.
        pgsacw(i) = min(psacws(i), cons*dtime*n0s(i)*qcic(i)*qcic(i)* &
                                   asn(i)*asn(i)/ (rho(i)*lams(i)**(bs+bs+two)))

!           if (pgsacw(i).lt.0_r8) then
!              write(iulog,*) "pgsacw,i,lams,cons",i,pgsacw(i),lams(i),cons
!           end if

! Mix rat converted into graupel as embryo (Reisner et al. 1998, orig M1990)
        dum= max(rhosn/(rhog-rhosn)*pgsacw(i), zero)

! Number concentraiton of embryo graupel from riming of snow
        nscng(i) = dum/mg0*rho(i)
! Limit max number converted to snow number  (dtime here correct? We think yes)
        nscng(i) = min(nscng(i),nsic(i)/dtime)

! Portion of riming left for snow
        psacws(i) = psacws(i) - pgsacw(i)
     else
        pgsacw(i) = zero
        nscng(i)  = zero
     end if

  enddo

end subroutine graupel_riming_liquid_snow

!========================================================================
!>\ingroup micro_mg_utils_mod
!!CHANGE IN Q,N COLLECTION RAIN BY GRAUPEL
subroutine graupel_collecting_rain(qric,qgic,umg,umr,ung,unr,rho,n0r,lamr,n0g,lamg,&
                                   pracg,npracg,mgncol)

  integer, intent(in) :: mgncol

  !MMR
  real(r8), dimension(mgncol), intent(in) :: qric  !rain mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qgic  !graupel mixing ratio

  !Mass weighted Fall speeds
  real(r8), dimension(mgncol), intent(in) :: umg ! graupel fall speed
  real(r8), dimension(mgncol), intent(in) :: umr ! rain fall speed

  !Number weighted fall speeds
  real(r8), dimension(mgncol), intent(in) :: ung ! graupel fall speed
  real(r8), dimension(mgncol), intent(in) :: unr ! rain fall speed

  real(r8), dimension(mgncol), intent(in) :: rho   ! air density

 ! Size parameters for rain
  real(r8), dimension(mgncol), intent(in) :: n0r
  real(r8), dimension(mgncol), intent(in) :: lamr

 ! Size parameters for graupel
  real(r8), dimension(mgncol), intent(in) :: n0g
  real(r8), dimension(mgncol), intent(in) :: lamg


  !Output process rates
  real(r8), dimension(mgncol), intent(out) :: pracg   ! Q collection rain by graupel
  real(r8), dimension(mgncol), intent(out) :: npracg  ! N collection rain by graupel

! Add collection of graupel by rain above freezing
! assume all rain collection by graupel above freezing is shed
! assume shed drops are 1 mm in size

  integer  :: i
  real(r8) :: cons41
  real(r8) :: cons32
  real(r8) :: dum, tx1, tx2, tx3, tx4, tx5, tx6

  cons41 = pi*pi*ecr*rhow
  cons32 = 0.5*pi*ecr

  do i=1,mgncol

     if (qric(i) >= 1.e-8_r8 .and. qgic(i) >= 1.e-8_r8) then

! pracg is mixing ratio of rain per sec collected by graupel/hail
        tx1 = 1.2_r8*umr(i) - 0.95_r8*umg(i)
        tx1 = sqrt(tx1*tx1+0.08_r8*umg(i)*umr(i))
        tx2 = 1.0_r8 / lamr(i)
        tx3 = 1.0_r8 / lamg(i)
        tx4 = tx2 * tx2
        tx5 = tx4 * tx4 * tx3
        tx6 = rho(i) * n0r(i) * n0g(i)


        pracg(i) = cons41 * tx1 * tx6 * tx5 * (5.0*tx4+tx3*(tx2+tx2+0.5*tx3))


!       pracg(i) = cons41*(((1.2_r8*umr(i)-0.95_r8*umg(i))**2._r8+      &
!            0.08_r8*umg(i)*umr(i))**0.5_r8*rho(i)*                     &
!            n0r(i)*n0g(i)/lamr(i)**3._r8*                              &
!            (5._r8/(lamr(i)**3._r8*lamg(i))+                           &
!            2._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+                     &
!            0.5_r8/(lamr(i)*lamg(i)**3._r8)))

! assume 1 mm drops are shed, get number shed per sec

        dum = pracg(i) / 5.2e-7_r8

        tx1 = unr(i) - ung(i)
        tx1 = sqrt(1.7_r8 * tx1 * tx1 + 0.3_r8*unr(i)*ung(i))
        tx4 = tx2 * tx3

        npracg(i) = cons32 * tx1 * tx6 * tx4 * (tx2*(tx2+tx3)+tx3*tx3)

!       npracg(i) = cons32*rho(i)*(1.7_r8*(unr(i)-ung(i))**2._r8+      &
!            0.3_r8*unr(i)*ung(i))**0.5_r8*n0r(i)*n0g(i)*              &
!            (1._r8/(lamr(i)**3._r8*lamg(i))+                          &
!             1._r8/(lamr(i)**2._r8*lamg(i)**2._r8)+                   &
!             1._r8/(lamr(i)*lamg(i)**3._r8))

! hm 7/15/13, remove limit so that the number of collected drops can smaller than
! number of shed drops
!            NPRACG(K)=MAX(NPRACG(K)-DUM,0.)
        npracg(i) = npracg(i) - dum
     else
        npracg(i) = zero
        pracg(i)  = zero
     end if

  enddo

end subroutine graupel_collecting_rain

!========================================================================
!>\ingroup micro_mg_utils_mod
!! Rain riming snow to graupel
!========================================================================
! Conversion of rimed rainwater onto snow converted to graupel
subroutine graupel_rain_riming_snow(pracs,npracs,psacr,qsic,qric,nric,nsic,n0s, &
                                    lams,n0r,lamr,dtime,pgracs,ngracs,mgncol)

  integer, intent(in) :: mgncol

  ! Accretion of rain by snow
  real(r8), dimension(mgncol), intent(inout) :: pracs
  real(r8), dimension(mgncol), intent(inout) :: npracs
  real(r8), dimension(mgncol), intent(inout) :: psacr  ! conversion due to coll of snow by rain


  real(r8), dimension(mgncol), intent(in) :: qsic  !snow mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qric  !rain mixing ratio

  real(r8), dimension(mgncol), intent(in) :: nric ! rain number concentration
  real(r8), dimension(mgncol), intent(in) :: nsic ! snow number concentration

 ! Size parameters for snow
  real(r8), dimension(mgncol), intent(in) :: n0s
  real(r8), dimension(mgncol), intent(in) :: lams

 ! Size parameters for rain
  real(r8), dimension(mgncol), intent(in) :: n0r
  real(r8), dimension(mgncol), intent(in) :: lamr

  real(r8),                    intent(in) :: dtime

  !Output process rates
  real(r8), dimension(mgncol), intent(out) :: pgracs  ! Q graupel due to collection rain by snow
  real(r8), dimension(mgncol), intent(out) :: ngracs  ! N graupel due to collection rain by snow

!Input: PRACS,NPRACS,PSACR,qni,qr,lams,lamr,nr,ns  Note: No PSACR in MG2
!Output:PGRACS,NGRACS,PRACS,PSACR

  integer  :: i
  real(r8) :: cons18
  real(r8) :: cons19
  real(r8) :: dum,fmult,tx1,tx2

  cons18 = rhosn*rhosn
  cons19 = rhow*rhow

  do i=1,mgncol

     fmult = zero

     if (pracs(i) > zero .and. qsic(i) >= 0.1e-3_r8 .and. qric(i) >= 0.1e-3_r8) then
        ! only allow conversion if qs > 0.1 and qr > 0.1 g/kg following rutledge and hobbs (1984)
        !if (qsic(i).ge.0.1e-3_r8.and.qric(i).ge.0.1e-3_r8) then
        ! portion of collected rainwater converted to graupel (reisner et al. 1998)
          tx1 = four / lams(i)
          tx2 = four / lamr(i)
          tx1 = tx1 * tx1 * tx1
          tx2 = tx2 * tx2 * tx2
          dum = cons18 * tx1 * tx1
          dum = one - max(zero, min(one, dum / (dum + cons19 * tx2 * tx2)))

!         dum = cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3 &
!              /(cons18*(4._r8/lams(i))**3*(4._r8/lams(i))**3+ &
!              cons19*(4._r8/lamr(i))**3*(4._r8/lamr(i))**3)
!         dum = min(dum,one)
!         dum = max(dum, zero)
!
!         pgracs(i) = (one-dum) * pracs(i)
!         ngracs(i) = (one-dum) * npracs(i)

          pgracs(i) = dum * pracs(i)
          ngracs(i) = dum * npracs(i)
          ! limit max number converted to min of either rain or snow number concentration
          ngracs(i) = min(ngracs(i),nric(i)/dtime)
          ngracs(i) = min(ngracs(i),nsic(i)/dtime)

          ! amount left for snow production
          pracs(i)  = pracs(i)  - pgracs(i)
          npracs(i) = npracs(i) - ngracs(i)

          ! conversion to graupel due to collection of snow by rain
!         psacr(i) = psacr(i) * (one-dum)
          psacr(i) = psacr(i) * dum
       else
          pgracs(i) = zero
          ngracs(i) = zero
        end if
     enddo

end subroutine graupel_rain_riming_snow

!========================================================================
! Rime Splintering
!========================================================================
!>\ingroup micro_mg_utils_mod
!! Rime splintering
subroutine graupel_rime_splintering(t,qcic,qric,qgic,psacwg,pracg,&
                                    qmultg,nmultg,qmultrg,nmultrg,mgncol)

  integer, intent(in) :: mgncol

  real(r8), dimension(mgncol), intent(in) :: t  !temperature

  !MMR
  real(r8), dimension(mgncol), intent(in) :: qcic  !liquid mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qric  !rain mixing ratio
  real(r8), dimension(mgncol), intent(in) :: qgic  !graupel mixing ratio

  ! Already calculated terms for collection
  real(r8), dimension(mgncol), intent(inout) :: psacwg ! collection droplets by graupel
  real(r8), dimension(mgncol), intent(inout) :: pracg  ! collection rain by graupel

  !Output process rates for splintering
  real(r8), dimension(mgncol), intent(out) :: qmultg  ! Q ice mult of droplets/graupel
  real(r8), dimension(mgncol), intent(out) :: nmultg  ! N ice mult of droplets/graupel
  real(r8), dimension(mgncol), intent(out) :: qmultrg  ! Q ice mult of rain/graupel
  real(r8), dimension(mgncol), intent(out) :: nmultrg  ! N ice mult of rain/graupel


!Input: qg,qc,qr, PSACWG,PRACG,T
!Output: NMULTG,QMULTG,NMULTRG,QMULTRG,PSACWG,PRACG

  integer  :: i
  real(r8) :: fmult
  real(r8) :: tm_3,tm_5,tm_8

  tm_3 = tmelt - 3._r8
  tm_5 = tmelt - 5._r8
  tm_8 = tmelt - 8._r8


!nmultg,qmultg                                                                             .
!========================================================================
  do i=1,mgncol

     nmultrg(i) = zero
     qmultrg(i) = zero
     nmultg(i)  = zero
     qmultg(i)  = zero

     if (qgic(i) >= 0.1e-3_r8) then
        if (qcic(i) >= 0.5e-3_r8 .or. qric(i) >= 0.1e-3_r8) then
           if (psacwg(i) > zero .or. pracg(i) > zero) then
              if (t(i) < tm_3 .and. t(i) > tm_8) then
                 if (t(i) > tm_3) then
                    fmult = zero
                 else if (t(i) <= tm_3 .and. t(i) > tm_5)  then
                    fmult = (tm_3-t(i)) * 0.5
                 else if (t(i) >= tm_8 .and. t(i) <= tm_5)   then
                    fmult = (t(i)-tm_8) * (one/three)
                 else if (t(i) < tm_8) then
                    fmult = zero
                 end if

! 1000 is to convert from kg to g  (Check if needed for MG)

! splintering from droplets accreted onto graupel

                 if (psacwg(i) > zero) then
                    nmultg(i) = 35.e4_r8*psacwg(i)*fmult*1000._r8
                    qmultg(i) = nmultg(i)*mmult

! constrain so that transfer of mass from graupel to ice cannot be more mass
! than was rimed onto graupel

                    qmultg(i) = min(qmultg(i),psacwg(i))
                    psacwg(i) = psacwg(i) - qmultg(i)

                 end if


!nmultrg,qmultrg                                                                             .
!========================================================================

! riming and splintering from accreted raindrops

! Factor of 1000. again (Check)

                 if (pracg(i) > zero) then
                    nmultrg(i) = 35.e4*pracg(i)*fmult*1000._r8
                    qmultrg(i) = nmultrg(i)*mmult

! constrain so that transfer of mass from graupel to ice cannot be more mass
! than was rimed onto graupel

                    qmultrg(i) = min(qmultrg(i),pracg(i))
                    pracg(i)   = pracg(i) - qmultrg(i)

                 end if

              end if
           end if
        end if
     end if
  enddo

end subroutine graupel_rime_splintering

!========================================================================
! Evaporation and sublimation of graupel
!========================================================================

!MERGE WITH RAIN AND SNOW EVAP
!
!subroutine graupel_sublimate_evap(t,q,qgic,rho,n0g,lamg,qvi,dv,mu,sc,bg,agn,&
!     prdg,eprdg,mgncol)
!
!  integer, intent(in) :: mgncol
!
!  real(r8), dimension(mgncol), intent(in) :: t  !temperature
!  real(r8), dimension(mgncol), intent(in) :: q  !specific humidity (mixing ratio)
!
!  !MMR
!  real(r8), dimension(mgncol), intent(in) :: qgic  !graupel mixing ratio
!
!  real(r8), dimension(mgncol), intent(in) :: rho   ! air density
!
! ! Size parameters for graupel
!  real(r8), dimension(mgncol), intent(in) :: n0g
!  real(r8), dimension(mgncol), intent(in) :: lamg
!
!  real(r8), dimension(mgncol), intent(in) :: qvi  !saturation humidity (ice)
!
!  real(r8), dimension(mgncol), intent(in) :: dv   ! water vapor diffusivity
!  real(r8), dimension(mgncol), intent(in) :: mu   ! viscosity
!  real(r8), dimension(mgncol), intent(in) :: sc   ! schmidt number
!
!  ! fallspeed parameters for graupel
!  ! Set AGN and BG  as input (in micro_mg3_0.F90)  (select hail or graupel as appropriate)
!  real(r8),                    intent(in) :: bg
!  real(r8), dimension(mgncol), intent(in) :: agn
!
!  ! Output tendencies (sublimation or evaporation of graupel)
!  real(r8), dimension(mgncol), intent(out) :: prdg
!  real(r8), dimension(mgncol), intent(out) :: eprdg
!
!  real(r8) :: cons11,cons36
!  real(r8) :: abi
!  real(r8) :: epsg
!  integer :: i
!
!  cons11=gamma(2.5_r8+bg/2._r8)  !bg will be different for graupel (bg) or hail (bh)
!  cons36=(2.5_r8+bg/2._r8)
!
!
!  do i=1,mgncol
!
!     abi = calc_ab(t(i), qvi(i), xxls)
!
!      if (qgic(i).ge.qsmall) then
!         epsg = 2._r8*pi*n0g(i)*rho(i)*dv(i)*                                &
!              (f1s/(lamg(i)*lamg(i))+                               &
!              f2s*(agn(i)*rho(i)/mu(i))**0.5_r8*                      &
!              sc(i)**(1._r8/3._r8)*cons11/                   &
!              (lamg(i)**cons36))
!      else
!         epsg = 0.
!      end if
!
!! vapor dpeosition on graupel
!      prdg(i) = epsg*(q(i)-qvi(i))/abi
!
!! make sure not pushed into ice supersat/subsat
!! put this in main mg3 code ... check for it ...
!! formula from reisner 2 scheme

!!
!!          dum = (qv3d(k)-qvi(k))/dt
!!
!!           fudgef = 0.9999
!!           sum_dep = prd(k)+prds(k)+mnuccd(k)+prdg(k)
!!
!!           if( (dum.gt.0. .and. sum_dep.gt.dum*fudgef) .or.                      &
!!               (dum.lt.0. .and. sum_dep.lt.dum*fudgef) ) then
!!             prdg(k) = fudgef*prdg(k)*dum/sum_dep
!!           endif
!
!! if cloud ice/snow/graupel vap deposition is neg, then assign to sublimation processes
!
!      eprdg(i)=0._r8
!
!      if (prdg(i).lt.0._r8) then
!         eprdg(i)=prdg(i)
!         prdg(i)=0.
!      end if
!
!   enddo
!
!end subroutine graupel_sublimate_evap

!========================================================================
!UTILITIES
!========================================================================

!>\ingroup micro_mg_utils_mod
pure function no_limiter()
  real(r8) :: no_limiter

  no_limiter = transfer(limiter_off, no_limiter)

end function no_limiter

!>\ingroup micro_mg_utils_mod
pure function limiter_is_on(lim)
  real(r8), intent(in) :: lim
  logical :: limiter_is_on

  limiter_is_on = transfer(lim, limiter_off) /= limiter_off

end function limiter_is_on

!>\ingroup micro_mg_utils_mod
FUNCTION gamma_incomp(muice, x)

  real(r8) :: gamma_incomp
  REAL(r8), intent(in) :: muice, x
  REAL(r8) :: xog, kg, alfa, auxx
  alfa = min(max(muice+1._r8, 1._r8), 20._r8)

  xog  = log(alfa -0.3068_r8)
  kg   = 1.44818_r8*(alfa**0.5357_r8)
  auxx = max(min(kg*(log(x)-xog), 30._r8), -30._r8)
  gamma_incomp = max(one/(one+exp(-auxx)), 1.0e-20) 

END FUNCTION gamma_incomp


end module micro_mg_utils
